diff --git a/packages/nox/src-thyra/NOX_Thyra_Group.C b/packages/nox/src-thyra/NOX_Thyra_Group.C index 73f55c19113b..8c20c9583e0d 100644 --- a/packages/nox/src-thyra/NOX_Thyra_Group.C +++ b/packages/nox/src-thyra/NOX_Thyra_Group.C @@ -380,7 +380,6 @@ NOX::Thyra::Group::getJacobian() const Teuchos::RCP< ::Thyra::PreconditionerBase > NOX::Thyra::Group::getNonconstPreconditioner() { - NOX_ASSERT(nonnull(prec_)); return prec_; } @@ -388,7 +387,6 @@ NOX::Thyra::Group::getNonconstPreconditioner() Teuchos::RCP > NOX::Thyra::Group::getPreconditioner() const { - NOX_ASSERT(nonnull(prec_)); return prec_; } diff --git a/packages/tempus/src/CMakeLists.txt b/packages/tempus/src/CMakeLists.txt index 9c34cab66c7b..5fd104fec5e1 100644 --- a/packages/tempus/src/CMakeLists.txt +++ b/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} ) + diff --git a/packages/tempus/src/Tempus_CombinedForwardSensitivityModelEvaluator.cpp b/packages/tempus/src/Tempus_CombinedForwardSensitivityModelEvaluator.cpp new file mode 100644 index 000000000000..62a1d99face4 --- /dev/null +++ b/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 diff --git a/packages/tempus/src/Tempus_CombinedForwardSensitivityModelEvaluator_decl.hpp b/packages/tempus/src/Tempus_CombinedForwardSensitivityModelEvaluator_decl.hpp new file mode 100644 index 000000000000..cdfa66a8445c --- /dev/null +++ b/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 +class CombinedForwardSensitivityModelEvaluator : + public Thyra::StateFuncModelEvaluatorBase { +public: + typedef Thyra::VectorBase Vector; + typedef Thyra::MultiVectorBase MultiVector; + + //! Constructor + /*! + * The optionally supplied parameter list supports the following options: + *
    + *
  • "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(). + *
  • "Sensitivity Parameter Index" (default: 0) Model evaluator + * parameter index for which sensitivities will be computed. + *
  • "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. + *
  • "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. + *
  • "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). + *
+ */ + CombinedForwardSensitivityModelEvaluator( + const Teuchos::RCP > & model, + const Teuchos::RCP& pList = Teuchos::null, + const Teuchos::RCP& dxdp_init = Teuchos::null, + const Teuchos::RCP& dx_dotdp_init = Teuchos::null, + const Teuchos::RCP& dx_dotdot_dp_init = Teuchos::null); + + //! Get the underlying model 'f' + Teuchos::RCP > getModel() const + { return model_; } + + /** \name Public functions overridden from ModelEvaulator. */ + //@{ + + Teuchos::RCP > get_p_space(int p) const; + + Teuchos::RCP > get_p_names(int p) const; + + Teuchos::RCP > get_x_space() const; + + Teuchos::RCP > get_f_space() const; + + Teuchos::RCP > create_W_op() const; + + Teuchos::RCP > + get_W_factory() const; + + Thyra::ModelEvaluatorBase::InArgs createInArgs() const; + + Thyra::ModelEvaluatorBase::InArgs getNominalValues() const; + + //@} + + static Teuchos::RCP getValidParameters(); + +private: + + Thyra::ModelEvaluatorBase::OutArgs createOutArgsImpl() const; + + void evalModelImpl( + const Thyra::ModelEvaluatorBase::InArgs &inArgs, + const Thyra::ModelEvaluatorBase::OutArgs &outArgs) const; + + + Thyra::ModelEvaluatorBase::InArgs prototypeInArgs_; + Thyra::ModelEvaluatorBase::OutArgs prototypeOutArgs_; + + Teuchos::RCP > model_; + Teuchos::RCP dxdp_init_; + Teuchos::RCP dx_dotdp_init_; + Teuchos::RCP 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 > dxdp_space_; + Teuchos::RCP > x_dxdp_space_; + Teuchos::RCP > dfdp_space_; + Teuchos::RCP > f_dfdp_space_; + + mutable Teuchos::RCP > my_dfdx_; + mutable Teuchos::RCP > my_dfdxdot_; + mutable Teuchos::RCP > my_dfdxdotdot_; +}; + +} // namespace Tempus + +#endif diff --git a/packages/tempus/src/Tempus_CombinedForwardSensitivityModelEvaluator_impl.hpp b/packages/tempus/src/Tempus_CombinedForwardSensitivityModelEvaluator_impl.hpp new file mode 100644 index 000000000000..b122f11644c2 --- /dev/null +++ b/packages/tempus/src/Tempus_CombinedForwardSensitivityModelEvaluator_impl.hpp @@ -0,0 +1,413 @@ +// @HEADER +// **************************************************************************** +// Tempus: Copyright (2017) Sandia Corporation +// +// Distributed under BSD 3-clause license (See accompanying file Copyright.txt) +// **************************************************************************** +// @HEADER + +#ifndef Tempus_CombinedForwardSensitivityModelEvaluator_impl_hpp +#define Tempus_CombinedForwardSensitivityModelEvaluator_impl_hpp + +#include "Thyra_DefaultMultiVectorProductVector.hpp" +#include "Thyra_MultiVectorLinearOp.hpp" +#include "Thyra_MultiVectorLinearOpWithSolveFactory.hpp" +#include "Thyra_ReuseLinearOpWithSolveFactory.hpp" + +namespace Tempus { + +template +CombinedForwardSensitivityModelEvaluator:: +CombinedForwardSensitivityModelEvaluator( + const Teuchos::RCP > & model, + const Teuchos::RCP& pList, + const Teuchos::RCP& dxdp_init, + const Teuchos::RCP& dx_dotdp_init, + const Teuchos::RCP& dx_dotdotdp_init) : + model_(model), + dxdp_init_(dxdp_init), + dx_dotdp_init_(dx_dotdp_init), + dx_dotdotdp_init_(dx_dotdotdp_init), + p_index_(0), + x_tangent_index_(1), + xdot_tangent_index_(2), + xdotdot_tangent_index_(3), + use_dfdp_as_tangent_(false) +{ + typedef Thyra::ModelEvaluatorBase MEB; + + // Set parameters + Teuchos::RCP pl = + Teuchos::rcp(new Teuchos::ParameterList); + if (pList != Teuchos::null) + *pl = *pList; + pl->validateParametersAndSetDefaults(*this->getValidParameters()); + use_dfdp_as_tangent_ = pl->get("Use DfDp as Tangent"); + p_index_ = pl->get("Sensitivity Parameter Index"); + x_tangent_index_ = pl->get("Sensitivity X Tangent Index"); + xdot_tangent_index_ = pl->get("Sensitivity X-Dot Tangent Index"); + xdotdot_tangent_index_ = pl->get("Sensitivity X-Dot-Dot Tangent Index"); + + num_param_ = model_->get_p_space(p_index_)->dim(); + dxdp_space_ = + Thyra::multiVectorProductVectorSpace(model_->get_x_space(), num_param_); + x_dxdp_space_ = + Thyra::multiVectorProductVectorSpace(model_->get_x_space(), num_param_+1); + dfdp_space_ = + Thyra::multiVectorProductVectorSpace(model_->get_f_space(), num_param_); + f_dfdp_space_ = + Thyra::multiVectorProductVectorSpace(model_->get_f_space(), num_param_+1); + + MEB::InArgs me_inArgs = model_->createInArgs(); + MEB::InArgsSetup inArgs; + inArgs.setModelEvalDescription(this->description()); + inArgs.setSupports(MEB::IN_ARG_x); + if (me_inArgs.supports(MEB::IN_ARG_x_dot)) + inArgs.setSupports(MEB::IN_ARG_x_dot); + if (me_inArgs.supports(MEB::IN_ARG_t)) + inArgs.setSupports(MEB::IN_ARG_t); + if (me_inArgs.supports(MEB::IN_ARG_alpha)) + inArgs.setSupports(MEB::IN_ARG_alpha); + if (me_inArgs.supports(MEB::IN_ARG_beta)) + inArgs.setSupports(MEB::IN_ARG_beta); + if (me_inArgs.supports(MEB::IN_ARG_W_x_dot_dot_coeff)) + inArgs.setSupports(MEB::IN_ARG_W_x_dot_dot_coeff); + + // Support additional parameters for x and xdot + inArgs.set_Np(me_inArgs.Np()); + prototypeInArgs_ = inArgs; + + MEB::OutArgs me_outArgs = model_->createOutArgs(); + MEB::OutArgsSetup outArgs; + outArgs.setModelEvalDescription(this->description()); + outArgs.set_Np_Ng(me_inArgs.Np(),0); + outArgs.setSupports(MEB::OUT_ARG_f); + if (me_outArgs.supports(MEB::OUT_ARG_W_op)) + outArgs.setSupports(MEB::OUT_ARG_W_op); + prototypeOutArgs_ = outArgs; + + TEUCHOS_ASSERT(me_outArgs.supports(MEB::OUT_ARG_DfDp, p_index_).supports(MEB::DERIV_MV_JACOBIAN_FORM)); + if (!use_dfdp_as_tangent_) + TEUCHOS_ASSERT(me_outArgs.supports(MEB::OUT_ARG_W_op)); +} + +template +Teuchos::RCP > +CombinedForwardSensitivityModelEvaluator:: +get_p_space(int p) const +{ + TEUCHOS_ASSERT(p < model_->Np()); + return model_->get_p_space(p); +} + +template +Teuchos::RCP > +CombinedForwardSensitivityModelEvaluator:: +get_p_names(int p) const +{ + TEUCHOS_ASSERT(p < model_->Np()); + return model_->get_p_names(p); +} + +template +Teuchos::RCP > +CombinedForwardSensitivityModelEvaluator:: +get_x_space() const +{ + return x_dxdp_space_; +} + +template +Teuchos::RCP > +CombinedForwardSensitivityModelEvaluator:: +get_f_space() const +{ + return f_dfdp_space_; +} + +template +Teuchos::RCP > +CombinedForwardSensitivityModelEvaluator:: +create_W_op() const +{ + Teuchos::RCP > op = model_->create_W_op(); + return Thyra::nonconstMultiVectorLinearOp(op, num_param_+1); +} + +template +Teuchos::RCP > +CombinedForwardSensitivityModelEvaluator:: +get_W_factory() const +{ + typedef Thyra::DefaultMultiVectorProductVectorSpace DMVPVS; + using Teuchos::RCP; + using Teuchos::rcp_dynamic_cast; + + Teuchos::RCP > factory = + model_->get_W_factory(); + RCP > mv_op = this->create_W_op(); + RCP mv_domain = + rcp_dynamic_cast(mv_op->domain()); + RCP mv_range = + rcp_dynamic_cast(mv_op->range()); + return Thyra::multiVectorLinearOpWithSolveFactory(factory, + mv_range, + mv_domain); +} + +template +Thyra::ModelEvaluatorBase::InArgs +CombinedForwardSensitivityModelEvaluator:: +createInArgs() const +{ + return prototypeInArgs_; +} + +template +Thyra::ModelEvaluatorBase::InArgs +CombinedForwardSensitivityModelEvaluator:: +getNominalValues() const +{ + typedef Thyra::ModelEvaluatorBase MEB; + typedef Thyra::DefaultMultiVectorProductVector DMVPV; + using Teuchos::RCP; + using Teuchos::rcp_dynamic_cast; + using Teuchos::Range1D; + + MEB::InArgs me_nominal = model_->getNominalValues(); + MEB::InArgs nominal = this->createInArgs(); + + const Teuchos::Range1D rng(1,num_param_); + const Scalar zero = Teuchos::ScalarTraits::zero(); + + // Set initial x. If dxdp_init == null, set initial dx/dp = 0 + RCP< const Thyra::VectorBase > me_x = me_nominal.get_x(); + if (me_x != Teuchos::null) { + RCP< Thyra::VectorBase > x = Thyra::createMember(*x_dxdp_space_); + RCP x_dxdp = rcp_dynamic_cast(x); + Thyra::assign(x_dxdp->getNonconstMultiVector()->col(0).ptr(), *me_x); + if (dxdp_init_ == Teuchos::null) + Thyra::assign(x_dxdp->getNonconstMultiVector()->subView(rng).ptr(), + zero); + else + Thyra::assign(x_dxdp->getNonconstMultiVector()->subView(rng).ptr(), + *dxdp_init_); + nominal.set_x(x); + } + + // Set initial xdot. If dx_dotdp_init == null, set initial dxdot/dp = 0 + RCP< const Thyra::VectorBase > me_xdot; + if (me_nominal.supports(MEB::IN_ARG_x_dot)) + me_xdot = me_nominal.get_x_dot(); + if (me_xdot != Teuchos::null) { + RCP< Thyra::VectorBase > xdot = Thyra::createMember(*x_dxdp_space_); + RCP xdot_dxdp = rcp_dynamic_cast(xdot); + Thyra::assign(xdot_dxdp->getNonconstMultiVector()->col(0).ptr(), *me_xdot); + if (dx_dotdp_init_ == Teuchos::null) + Thyra::assign(xdot_dxdp->getNonconstMultiVector()->subView(rng).ptr(), + zero); + else + Thyra::assign(xdot_dxdp->getNonconstMultiVector()->subView(rng).ptr(), + *dx_dotdp_init_); + nominal.set_x_dot(xdot); + } + + // Set initial xdotdot. If dx_dotdotdp_init == null, set initial dxdotdot/dp = 0 + RCP< const Thyra::VectorBase > me_xdotdot; + if (me_nominal.supports(MEB::IN_ARG_x_dot_dot)) + me_xdotdot = me_nominal.get_x_dot_dot(); + if (me_xdotdot != Teuchos::null) { + RCP< Thyra::VectorBase > xdotdot = + Thyra::createMember(*x_dxdp_space_); + RCP xdotdot_dxdp = rcp_dynamic_cast(xdotdot); + Thyra::assign(xdotdot_dxdp->getNonconstMultiVector()->col(0).ptr(), *me_xdotdot); + if (dx_dotdotdp_init_ == Teuchos::null) + Thyra::assign(xdotdot_dxdp->getNonconstMultiVector()->subView(rng).ptr(), + zero); + else + Thyra::assign(xdotdot_dxdp->getNonconstMultiVector()->subView(rng).ptr(), + *dx_dotdotdp_init_); + nominal.set_x_dot_dot(xdotdot); + } + + const int np = model_->Np(); + for (int i=0; i +Thyra::ModelEvaluatorBase::OutArgs +CombinedForwardSensitivityModelEvaluator:: +createOutArgsImpl() const +{ + return prototypeOutArgs_; +} + +template +void +CombinedForwardSensitivityModelEvaluator:: +evalModelImpl(const Thyra::ModelEvaluatorBase::InArgs &inArgs, + const Thyra::ModelEvaluatorBase::OutArgs &outArgs) const +{ + typedef Thyra::ModelEvaluatorBase MEB; + typedef Thyra::DefaultMultiVectorProductVector DMVPV; + using Teuchos::RCP; + using Teuchos::rcp_dynamic_cast; + using Teuchos::Range1D; + + // setup input arguments for model + RCP< const Thyra::VectorBase > x, xdot, xdotdot; + RCP< const Thyra::MultiVectorBase > dxdp, dxdotdp, dxdotdotdp; + MEB::InArgs me_inArgs = model_->getNominalValues(); + RCP x_dxdp = rcp_dynamic_cast(inArgs.get_x()); + x = x_dxdp->getMultiVector()->col(0); + dxdp = x_dxdp->getMultiVector()->subView(Range1D(1,num_param_)); + me_inArgs.set_x(x); + if (use_dfdp_as_tangent_) { + RCP< const Thyra::VectorBase > dxdp_vec = + Thyra::multiVectorProductVector(dxdp_space_, dxdp); + me_inArgs.set_p(x_tangent_index_, dxdp_vec); + } + if (me_inArgs.supports(MEB::IN_ARG_x_dot)) { + if (inArgs.get_x_dot() != Teuchos::null) { + RCP xdot_dxdotdp = + rcp_dynamic_cast(inArgs.get_x_dot()); + xdot = xdot_dxdotdp->getMultiVector()->col(0); + dxdotdp = xdot_dxdotdp->getMultiVector()->subView(Range1D(1,num_param_)); + me_inArgs.set_x_dot(xdot); + if (use_dfdp_as_tangent_) { + RCP< const Thyra::VectorBase > dxdotdp_vec = + Thyra::multiVectorProductVector(dxdp_space_, dxdotdp); + me_inArgs.set_p(xdot_tangent_index_, dxdotdp_vec); + } + } + else // clear out xdot if it was set in nominalValues + me_inArgs.set_x_dot(Teuchos::null); + } + if (me_inArgs.supports(MEB::IN_ARG_x_dot_dot)) { + if (inArgs.get_x_dot_dot() != Teuchos::null) { + RCP xdotdot_dxdotdotdp = + rcp_dynamic_cast(inArgs.get_x_dot_dot()); + xdotdot = xdotdot_dxdotdotdp->getMultiVector()->col(0); + dxdotdotdp = xdotdot_dxdotdotdp->getMultiVector()->subView(Range1D(1,num_param_)); + me_inArgs.set_x_dot_dot(xdotdot); + if (use_dfdp_as_tangent_) { + RCP< const Thyra::VectorBase > dxdotdotdp_vec = + Thyra::multiVectorProductVector(dxdp_space_, dxdotdotdp); + me_inArgs.set_p(xdotdot_tangent_index_, dxdotdotdp_vec); + } + } + else // clear out xdotdot if it was set in nominalValues + me_inArgs.set_x_dot_dot(Teuchos::null); + } + if (me_inArgs.supports(MEB::IN_ARG_t)) + me_inArgs.set_t(inArgs.get_t()); + if (me_inArgs.supports(MEB::IN_ARG_alpha)) + me_inArgs.set_alpha(inArgs.get_alpha()); + if (me_inArgs.supports(MEB::IN_ARG_beta)) + me_inArgs.set_beta(inArgs.get_beta()); + if (me_inArgs.supports(MEB::IN_ARG_W_x_dot_dot_coeff)) + me_inArgs.set_W_x_dot_dot_coeff(inArgs.get_W_x_dot_dot_coeff()); + + // Set parameters -- be careful to only set ones that were set in our + // inArgs to not null out any specified through nominalValues or + // dx/dp above + const int np = me_inArgs.Np(); + for (int i=0; i > f; + RCP< Thyra::MultiVectorBase > dfdp; + MEB::OutArgs me_outArgs = model_->createOutArgs(); + if (outArgs.get_f() != Teuchos::null) { + RCP f_dfdp = rcp_dynamic_cast(outArgs.get_f()); + f = f_dfdp->getNonconstMultiVector()->col(0); + dfdp = f_dfdp->getNonconstMultiVector()->subView(Range1D(1,num_param_)); + me_outArgs.set_f(f); + me_outArgs.set_DfDp(p_index_, dfdp); + } + if (outArgs.get_W_op() != Teuchos::null) { + RCP< Thyra::LinearOpBase > op = outArgs.get_W_op(); + RCP< Thyra::MultiVectorLinearOp > mv_op = + rcp_dynamic_cast< Thyra::MultiVectorLinearOp >(op); + me_outArgs.set_W_op(mv_op->getNonconstLinearOp()); + } + + // build residual and jacobian + model_->evalModel(me_inArgs, me_outArgs); + + // Compute (df/dx) * (dx/dp) + (df/dxdot) * (dxdot/dp) + (df/dxdotdot) * (dxdotdot/dp) + (df/dp) + // if the underlying ME doesn't already do this. This requires computing + // df/dx, df/dxdot, df/dxdotdot as separate operators + if (!use_dfdp_as_tangent_) { + if (dxdp != Teuchos::null && dfdp != Teuchos::null) { + if (my_dfdx_ == Teuchos::null) + my_dfdx_ = model_->create_W_op(); + MEB::OutArgs meo = model_->createOutArgs(); + meo.set_W_op(my_dfdx_); + if (me_inArgs.supports(MEB::IN_ARG_alpha)) + me_inArgs.set_alpha(0.0); + if (me_inArgs.supports(MEB::IN_ARG_beta)) + me_inArgs.set_beta(1.0); + if (me_inArgs.supports(MEB::IN_ARG_W_x_dot_dot_coeff)) + me_inArgs.set_W_x_dot_dot_coeff(0.0); + model_->evalModel(me_inArgs, meo); + my_dfdx_->apply(Thyra::NOTRANS, *dxdp, dfdp.ptr(), Scalar(1.0), Scalar(1.0)); + } + if (dxdotdp != Teuchos::null && dfdp != Teuchos::null) { + if (my_dfdxdot_ == Teuchos::null) + my_dfdxdot_ = model_->create_W_op(); + MEB::OutArgs meo = model_->createOutArgs(); + meo.set_W_op(my_dfdxdot_); + if (me_inArgs.supports(MEB::IN_ARG_alpha)) + me_inArgs.set_alpha(1.0); + if (me_inArgs.supports(MEB::IN_ARG_beta)) + me_inArgs.set_beta(0.0); + if (me_inArgs.supports(MEB::IN_ARG_W_x_dot_dot_coeff)) + me_inArgs.set_W_x_dot_dot_coeff(0.0); + model_->evalModel(me_inArgs, meo); + my_dfdxdot_->apply(Thyra::NOTRANS, *dxdotdp, dfdp.ptr(), Scalar(1.0), Scalar(1.0)); + } + if (dxdotdotdp != Teuchos::null && dfdp != Teuchos::null) { + if (my_dfdxdotdot_ == Teuchos::null) + my_dfdxdotdot_ = model_->create_W_op(); + MEB::OutArgs meo = model_->createOutArgs(); + meo.set_W_op(my_dfdxdotdot_); + if (me_inArgs.supports(MEB::IN_ARG_alpha)) + me_inArgs.set_alpha(0.0); + if (me_inArgs.supports(MEB::IN_ARG_beta)) + me_inArgs.set_beta(0.0); + if (me_inArgs.supports(MEB::IN_ARG_W_x_dot_dot_coeff)) + me_inArgs.set_W_x_dot_dot_coeff(1.0); + model_->evalModel(me_inArgs, meo); + my_dfdxdotdot_->apply(Thyra::NOTRANS, *dxdotdotdp, dfdp.ptr(), Scalar(1.0), Scalar(1.0)); + } + } +} + +template +Teuchos::RCP +CombinedForwardSensitivityModelEvaluator:: +getValidParameters() +{ + Teuchos::RCP pl = Teuchos::parameterList(); + pl->set("Use DfDp as Tangent", false); + pl->set("Sensitivity Parameter Index", 0); + pl->set("Sensitivity X Tangent Index", 1); + pl->set("Sensitivity X-Dot Tangent Index", 2); + pl->set("Sensitivity X-Dot-Dot Tangent Index", 3); + return pl; +} + +} // namespace Tempus + +#endif diff --git a/packages/tempus/src/Tempus_IntegratorBasic_decl.hpp b/packages/tempus/src/Tempus_IntegratorBasic_decl.hpp index 328dce3d5bf3..ebe1e5b89e81 100644 --- a/packages/tempus/src/Tempus_IntegratorBasic_decl.hpp +++ b/packages/tempus/src/Tempus_IntegratorBasic_decl.hpp @@ -100,9 +100,9 @@ class IntegratorBasic : virtual public Tempus::Integrator Teuchos::RCP > state = Teuchos::null); /// Set the initial state from Thyra::VectorBase(s) virtual void setInitialState(Scalar t0, - Teuchos::RCP > x0, - Teuchos::RCP > xdot0 = Teuchos::null, - Teuchos::RCP > xdotdot0 = Teuchos::null); + Teuchos::RCP > x0, + Teuchos::RCP > xdot0 = Teuchos::null, + Teuchos::RCP > xdotdot0 = Teuchos::null); /// Get the SolutionHistory virtual Teuchos::RCP > getSolutionHistory() const override {return solutionHistory_;} @@ -138,6 +138,9 @@ class IntegratorBasic : virtual public Tempus::Integrator /// Get current state virtual Teuchos::RCP > getCurrentState() {return solutionHistory_->getCurrentState();} + + Teuchos::RCP getIntegratorParameterList() + { return integratorPL_; } //@} /// Parse when screen output should be executed diff --git a/packages/tempus/src/Tempus_IntegratorBasic_impl.hpp b/packages/tempus/src/Tempus_IntegratorBasic_impl.hpp index 212569969793..6ca4c9be05c1 100644 --- a/packages/tempus/src/Tempus_IntegratorBasic_impl.hpp +++ b/packages/tempus/src/Tempus_IntegratorBasic_impl.hpp @@ -156,9 +156,9 @@ setInitialState(Teuchos::RCP > state) template void IntegratorBasic:: setInitialState(Scalar t0, - Teuchos::RCP > x0, - Teuchos::RCP > xdot0, - Teuchos::RCP > xdotdot0) + Teuchos::RCP > x0, + Teuchos::RCP > xdot0, + Teuchos::RCP > xdotdot0) { using Teuchos::RCP; using Teuchos::ParameterList; @@ -192,7 +192,7 @@ setInitialState(Scalar t0, Thyra::assign(xdotdot.ptr(), *(xdotdot0)); RCP > newState = rcp(new SolutionState( - md, x0, xdot, xdotdot, stepper_->getDefaultStepperState())); + md, x0->clone_v(), xdot, xdotdot, stepper_->getDefaultStepperState())); solutionHistory_->addState(newState); } @@ -699,6 +699,17 @@ Teuchos::RCP > integratorBasic( return(integrator); } +/// Non-member constructor +template +Teuchos::RCP > integratorBasic( + Teuchos::RCP pList, + const Teuchos::RCP >& stepper) +{ + Teuchos::RCP > integrator = + Teuchos::rcp(new Tempus::IntegratorBasic(pList, stepper)); + return(integrator); +} + /// Non-member constructor template Teuchos::RCP > integratorBasic( diff --git a/packages/tempus/src/Tempus_IntegratorForwardSensitivity.cpp b/packages/tempus/src/Tempus_IntegratorForwardSensitivity.cpp new file mode 100644 index 000000000000..aa06d090adbf --- /dev/null +++ b/packages/tempus/src/Tempus_IntegratorForwardSensitivity.cpp @@ -0,0 +1,37 @@ +// @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_IntegratorForwardSensitivity.hpp" +#include "Tempus_IntegratorForwardSensitivity_impl.hpp" + +namespace Tempus { + + TEMPUS_INSTANTIATE_TEMPLATE_CLASS(IntegratorForwardSensitivity) + + // non-member ctor + template Teuchos::RCP > + integratorForwardSensitivity( + Teuchos::RCP parameterList, + const Teuchos::RCP >& model); + + // non-member ctor + template Teuchos::RCP > + integratorForwardSensitivity( + const Teuchos::RCP >& model, + std::string stepperType); + + // non-member ctor + template Teuchos::RCP > + integratorForwardSensitivity(); + +} // namespace Tempus + +#endif diff --git a/packages/tempus/src/Tempus_IntegratorForwardSensitivity_decl.hpp b/packages/tempus/src/Tempus_IntegratorForwardSensitivity_decl.hpp new file mode 100644 index 000000000000..7ca8ef85e6ec --- /dev/null +++ b/packages/tempus/src/Tempus_IntegratorForwardSensitivity_decl.hpp @@ -0,0 +1,259 @@ +// @HEADER +// **************************************************************************** +// Tempus: Copyright (2017) Sandia Corporation +// +// Distributed under BSD 3-clause license (See accompanying file Copyright.txt) +// **************************************************************************** +// @HEADER + +#ifndef Tempus_IntegratorForwardSensitivity_decl_hpp +#define Tempus_IntegratorForwardSensitivity_decl_hpp + +// Tempus +#include "Tempus_IntegratorBasic.hpp" +#include "Tempus_StepperStaggeredForwardSensitivity.hpp" + +namespace Tempus { + + +/** \brief Time integrator implementing forward sensitivity analysis */ +/** + * This integrator implements forward parameter sensitivity analysis by + * propagating the derivative of the solution with respect to model parameters + * alongside the solution. It supports sensitivity propagation methods: + *
    + *
  • "Combined" where the sensitivity and state equations are solved + * simultaneously. This is most appropriate for explicit time + * integration methods or implicit methods where a very-lightweight + * nonlinear solution strategy is used. + *
  • "Staggered" where the sensitivity equations are solved at each time + * step immediately after the state equations are solved. This is + * useful for implicit methods since it saves the cost of solving the + * sensitivity equations while the state equations are being solved. It + * generally does not work for explicit methods (and wouldn't be more + * efficient than the combined method even if it did). + *
+ * + * Note that this integrator implements all of the same functions as the + * IntegratorBasic, but is not derived from IntegratorBasic. It also provides + * functions for setting the sensitivity initial conditions and extracting the + * sensitivity at the final time. Also the vectors stored in the solution + * history store product vectors of the state and sensitivities using + * Thyra;:DefaultMultiVectorProductVector. + */ +template +class IntegratorForwardSensitivity : virtual public Tempus::Integrator +{ +public: + + /** \brief Constructor with ParameterList and model, and will be fully + * initialized. */ + /*! + * In addition to all of the regular integrator options, the supplied + * parameter list supports the following options contained within a sublist + * "Sensitivities" from the top-level parameter list: + *
    + *
  • "Sensitivity Method" (default: "Combined") The sensitivity + * analysis method as described above. + *
  • "Reuse State Linear Solver" (default: false) For the staggered + * method, whether to reuse the model's W matrix, solver, and + * preconditioner when solving the sensitivity equations. If they + * can be reused, substantial savings in compute time are possible. + *
  • "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 in the Tempus::CombinedForwardSensitivityModelEvaluator + * documentation. + *
  • "Sensitivity Parameter Index" (default: 0) Model evaluator + * parameter index for which sensitivities will be computed. + *
  • "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. + *
  • "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. + *
  • "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). + *
+ */ + IntegratorForwardSensitivity( + Teuchos::RCP pList, + const Teuchos::RCP >& model); + + /** \brief Constructor with model and "Stepper Type" and is fully initialized with default settings. */ + IntegratorForwardSensitivity( + const Teuchos::RCP >& model, + std::string stepperType); + + /// Destructor + /** \brief Constructor that requires a subsequent setParameterList, setStepper, and initialize calls. */ + IntegratorForwardSensitivity(); + + /// Destructor + virtual ~IntegratorForwardSensitivity() {} + + /// \name Basic integrator methods + //@{ + /// Advance the solution to timeMax, and return true if successful. + virtual bool advanceTime() + { return integrator_->advanceTime(); } + /// Advance the solution to timeFinal, and return true if successful. + virtual bool advanceTime(const Scalar timeFinal) override + { return integrator_->advanceTime(timeFinal); } + /// Perform tasks before start of integrator. + virtual void startIntegrator() + { integrator_->startIntegrator(); } + /// Start time step. + virtual void startTimeStep() + { integrator_->startTimeStep(); } + /// Only accept step after meeting time step criteria. + virtual void acceptTimeStep() + { integrator_->acceptTimeStep(); } + /// Perform tasks after end of integrator. + virtual void endIntegrator() + { integrator_->endIntegrator(); } + /// Return a copy of the Tempus ParameterList + virtual Teuchos::RCP getTempusParameterList() override + { return integrator_->getTempusParameterList(); } + virtual void setTempusParameterList(Teuchos::RCP pl) override + { integrator_->setTempusParameterList(pl); } + //@} + + /// \name Accessor methods + //@{ + /// Get current time + virtual Scalar getTime() const override + { return integrator_->getTime(); } + /// Get current index + virtual Scalar getIndex() const override + { return integrator_->getIndex(); } + /// Get Status + virtual Status getStatus() const override + { return integrator_->getStatus(); } + /// Get the Stepper + virtual Teuchos::RCP > getStepper() const override + { return integrator_->getStepper(); } + + /// Set the Stepper + virtual void setStepper(Teuchos::RCP > model); + + /// Set the Stepper + virtual void setStepperWStepper(Teuchos::RCP > stepper) + { integrator_->setStepperWStepper(stepper); } + /// Set the initial state which has the initial conditions + virtual void setInitialState( + Teuchos::RCP > state = Teuchos::null) + { integrator_->setInitialState(state); } + + /// Set the initial state from Thyra::VectorBase(s) + virtual void setInitialState( + Scalar t0, + Teuchos::RCP > x0, + Teuchos::RCP > xdot0 = Teuchos::null, + Teuchos::RCP > xdotdot0 = Teuchos::null, + Teuchos::RCP > DxDp0 = Teuchos::null, + Teuchos::RCP > DxdotDp0 = Teuchos::null, + Teuchos::RCP > DxdotdotDp0 = Teuchos::null); + + /// Get the SolutionHistory + virtual Teuchos::RCP > getSolutionHistory() const override + { return integrator_->getSolutionHistory(); } + /// Set the SolutionHistory + virtual void setSolutionHistory( + Teuchos::RCP > sh = Teuchos::null) + { integrator_->setSolutionHistory(sh); } + /// Get the TimeStepControl + virtual Teuchos::RCP > getTimeStepControl() const override + { return integrator_->getTimeStepControl(); } + /// Set the TimeStepControl + virtual void setTimeStepControl( + Teuchos::RCP > tsc = Teuchos::null) + { integrator_->setTimeStepControl(tsc); } + /// Get the Observer + virtual Teuchos::RCP > getObserver() + { return integrator_->getObserver(); } + /// Set the Observer + virtual void setObserver( + Teuchos::RCP > obs = Teuchos::null) + { integrator_->setObserver(obs); } + /// Initializes the Integrator after set* function calls + virtual void initialize() + { integrator_->initialize(); } + + /// Get current the solution, x + virtual Teuchos::RCP > getX() const; + virtual Teuchos::RCP > getDxDp() const; + /// Get current the time derivative of the solution, xdot + virtual Teuchos::RCP > getXdot() const; + virtual Teuchos::RCP > getDxdotDp() const; + /// Get current the second time derivative of the solution, xdotdot + virtual Teuchos::RCP > getXdotdot() const; + virtual Teuchos::RCP > getDxdotdotDp() const; + + /// Get current state + virtual Teuchos::RCP > getCurrentState() + {return integrator_->getCurrentState();} + //@} + + /// Parse when screen output should be executed + void parseScreenOutput() { integrator_->parseScreenOutput(); } + + /// \name Overridden from Teuchos::ParameterListAcceptor + //@{ + void setParameterList(const Teuchos::RCP & pl); + Teuchos::RCP getNonconstParameterList() override + { return tempus_pl_; } + Teuchos::RCP unsetParameterList() override; + + Teuchos::RCP getValidParameters() + const override; + //@} + + /// \name Overridden from Teuchos::Describable + //@{ + std::string description() const override; + void describe(Teuchos::FancyOStream & out, + const Teuchos::EVerbosityLevel verbLevel) const override; + //@} + +protected: + + // Create sensitivity model evaluator from application model + void + createSensitivityModelAndStepper( + const Teuchos::RCP >& model); + + Teuchos::RCP > model_; + Teuchos::RCP > sens_model_; + Teuchos::RCP > sens_stepper_; + Teuchos::RCP > integrator_; + Teuchos::RCP tempus_pl_; + Teuchos::RCP sens_pl_; + Teuchos::RCP stepper_pl_; + bool use_combined_method_; +}; + +/// Non-member constructor +template +Teuchos::RCP > +integratorForwardSensitivity( + Teuchos::RCP pList, + const Teuchos::RCP >& model); + +/// Non-member constructor +template +Teuchos::RCP > +integratorForwardSensitivity( + const Teuchos::RCP >& model, + std::string stepperType); + +/// Non-member constructor +template +Teuchos::RCP > +integratorForwardSensitivity(); + +} // namespace Tempus + +#endif // Tempus_IntegratorForwardSensitivity_decl_hpp diff --git a/packages/tempus/src/Tempus_IntegratorForwardSensitivity_impl.hpp b/packages/tempus/src/Tempus_IntegratorForwardSensitivity_impl.hpp new file mode 100644 index 000000000000..16ec170458ab --- /dev/null +++ b/packages/tempus/src/Tempus_IntegratorForwardSensitivity_impl.hpp @@ -0,0 +1,359 @@ +// @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 "Teuchos_VerboseObjectParameterListHelpers.hpp" +#include "Thyra_DefaultMultiVectorProductVector.hpp" +#include "Tempus_CombinedForwardSensitivityModelEvaluator.hpp" + +namespace Tempus { + +template +IntegratorForwardSensitivity:: +IntegratorForwardSensitivity( + Teuchos::RCP inputPL, + const Teuchos::RCP >& model) +{ + model_ = model; + integrator_ = integratorBasic(); + this->setParameterList(inputPL); + createSensitivityModelAndStepper(model); + if (use_combined_method_) + integrator_ = integratorBasic(tempus_pl_, sens_model_); + else { + integrator_ = integratorBasic(); + integrator_->setTempusParameterList(tempus_pl_); + integrator_->setStepperWStepper(sens_stepper_); + integrator_->initialize(); + } +} + +template +IntegratorForwardSensitivity:: +IntegratorForwardSensitivity( + const Teuchos::RCP >& model, + std::string stepperType) +{ + model_ = model; + integrator_ = integratorBasic(); + this->setParameterList(Teuchos::null); + createSensitivityModelAndStepper(model); + if (use_combined_method_) + integrator_ = integratorBasic(sens_model_, stepperType); + else { + integrator_ = integratorBasic(); + integrator_->setParameterList(tempus_pl_); + integrator_->setStepperWStepper(sens_stepper_); + integrator_->initialize(); + } + +} + +template +IntegratorForwardSensitivity:: +IntegratorForwardSensitivity() +{ + integrator_ = integratorBasic(); + this->setParameterList(Teuchos::null); +} + +template +void IntegratorForwardSensitivity:: +setStepper( + Teuchos::RCP > model) +{ + createSensitivityModelAndStepper(model); + if (use_combined_method_) + integrator_->setStepper(sens_model_); + else + integrator_->setStepperWStepper(sens_stepper_); +} + +template +void IntegratorForwardSensitivity:: +setInitialState(Scalar t0, + Teuchos::RCP > x0, + Teuchos::RCP > xdot0, + Teuchos::RCP > xdotdot0, + Teuchos::RCP > DxDp0, + Teuchos::RCP > DxdotDp0, + Teuchos::RCP > DxdotdotDp0) +{ + using Teuchos::RCP; + using Teuchos::rcp_dynamic_cast; + using Thyra::VectorSpaceBase; + using Thyra::assign; + using Thyra::createMember; + typedef Thyra::DefaultMultiVectorProductVector DMVPV; + + // + // Create and initialize product X, Xdot, Xdotdot + + RCP< const VectorSpaceBase > space; + if (use_combined_method_) + space = sens_model_->get_x_space(); + else + space = sens_stepper_->get_x_space(); + RCP X = rcp_dynamic_cast(createMember(space)); + RCP Xdot = rcp_dynamic_cast(createMember(space)); + RCP Xdotdot = rcp_dynamic_cast(createMember(space)); + + const int num_param = X->getNonconstMultiVector()->domain()->dim()-1; + const Scalar zero = Teuchos::ScalarTraits::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_->setInitialState(t0, X, Xdot, Xdotdot); +} + +template +Teuchos::RCP > +IntegratorForwardSensitivity:: +getX() const +{ + using Teuchos::RCP; + using Teuchos::rcp_dynamic_cast; + typedef Thyra::DefaultMultiVectorProductVector DMVPV; + + RCP X = rcp_dynamic_cast(integrator_->getX()); + return X->getMultiVector()->col(0); +} + +template +Teuchos::RCP > +IntegratorForwardSensitivity:: +getDxDp() const +{ + using Teuchos::RCP; + using Teuchos::rcp_dynamic_cast; + typedef Thyra::DefaultMultiVectorProductVector DMVPV; + + RCP X = rcp_dynamic_cast(integrator_->getX()); + const int num_param = X->getMultiVector()->domain()->dim()-1; + const Teuchos::Range1D rng(1,num_param); + return X->getMultiVector()->subView(rng); +} + +template +Teuchos::RCP > +IntegratorForwardSensitivity:: +getXdot() const +{ + using Teuchos::RCP; + using Teuchos::rcp_dynamic_cast; + typedef Thyra::DefaultMultiVectorProductVector DMVPV; + + RCP Xdot = rcp_dynamic_cast(integrator_->getXdot()); + return Xdot->getMultiVector()->col(0); +} + +template +Teuchos::RCP > +IntegratorForwardSensitivity:: +getDxdotDp() const +{ + using Teuchos::RCP; + using Teuchos::rcp_dynamic_cast; + typedef Thyra::DefaultMultiVectorProductVector DMVPV; + + RCP Xdot = rcp_dynamic_cast(integrator_->getXdot()); + const int num_param = Xdot->getMultiVector()->domain()->dim()-1; + const Teuchos::Range1D rng(1,num_param); + return Xdot->getMultiVector()->subView(rng); +} + +template +Teuchos::RCP > +IntegratorForwardSensitivity:: +getXdotdot() const +{ + using Teuchos::RCP; + using Teuchos::rcp_dynamic_cast; + typedef Thyra::DefaultMultiVectorProductVector DMVPV; + + RCP Xdotdot = + rcp_dynamic_cast(integrator_->getXdotdot()); + return Xdotdot->getMultiVector()->col(0); +} + +template +Teuchos::RCP > +IntegratorForwardSensitivity:: +getDxdotdotDp() const +{ + using Teuchos::RCP; + using Teuchos::rcp_dynamic_cast; + typedef Thyra::DefaultMultiVectorProductVector DMVPV; + + RCP Xdotdot = + rcp_dynamic_cast(integrator_->getXdotdot()); + const int num_param = Xdotdot->getMultiVector()->domain()->dim()-1; + const Teuchos::Range1D rng(1,num_param); + return Xdotdot->getMultiVector()->subView(rng); +} + +template +std::string +IntegratorForwardSensitivity:: +description() const +{ + std::string name = "Tempus::IntegratorForwardSensitivity"; + return(name); +} + +template +void +IntegratorForwardSensitivity:: +describe( + Teuchos::FancyOStream &out, + const Teuchos::EVerbosityLevel verbLevel) const +{ + out << description() << "::describe" << std::endl; + integrator_->describe(out, verbLevel); +} + +template +void +IntegratorForwardSensitivity:: +setParameterList(const Teuchos::RCP & inputPL) +{ + tempus_pl_ = Teuchos::parameterList(); + if (inputPL != Teuchos::null) + *tempus_pl_ = *inputPL; + tempus_pl_->setParametersNotAlreadySet(*this->getValidParameters()); + sens_pl_ = Teuchos::sublist(tempus_pl_, "Sensitivities", false); + std::string integratorName = + tempus_pl_->get("Integrator Name", "Default Integrator"); + std::string stepperName = + tempus_pl_->sublist(integratorName).get("Stepper Name"); + stepper_pl_ = Teuchos::sublist(tempus_pl_, stepperName, true); + use_combined_method_ = + sens_pl_->get("Sensitivity Method") == "Combined"; +} + +template +Teuchos::RCP +IntegratorForwardSensitivity:: +unsetParameterList() +{ + Teuchos::RCP temp_param_list = tempus_pl_; + tempus_pl_ = Teuchos::null; + sens_pl_ = Teuchos::null; + stepper_pl_ = Teuchos::null; + integrator_->unsetParameterList(); + return temp_param_list; +} + +template +Teuchos::RCP +IntegratorForwardSensitivity:: +getValidParameters() const +{ + Teuchos::RCP pl = + Teuchos::rcp(new Teuchos::ParameterList); + Teuchos::RCP integrator_pl = + integrator_->getValidParameters(); + Teuchos::RCP sensitivity_pl = + CombinedForwardSensitivityModelEvaluator::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); + + return pl; +} + +template +void +IntegratorForwardSensitivity:: +createSensitivityModelAndStepper( + const Teuchos::RCP >& model) +{ + using Teuchos::rcp; + + Teuchos::RCP spl = Teuchos::parameterList(); + *spl = *sens_pl_; + spl->remove("Sensitivity Method"); + + if (use_combined_method_) { + spl->remove("Reuse State Linear Solver"); + sens_model_ = + rcp(new CombinedForwardSensitivityModelEvaluator(model_, spl)); + } + else { + sens_stepper_ = + rcp(new StepperStaggeredForwardSensitivity( + model_, stepper_pl_, spl)); + } +} + +/// Non-member constructor +template +Teuchos::RCP > +integratorForwardSensitivity( + Teuchos::RCP pList, + const Teuchos::RCP >& model) +{ + Teuchos::RCP > integrator = + Teuchos::rcp(new Tempus::IntegratorForwardSensitivity(pList, model)); + return(integrator); +} + +/// Non-member constructor +template +Teuchos::RCP > +integratorForwardSensitivity( + const Teuchos::RCP >& model, + std::string stepperType) +{ + Teuchos::RCP > integrator = + Teuchos::rcp(new Tempus::IntegratorForwardSensitivity(model, stepperType)); + return(integrator); +} + +/// Non-member constructor +template +Teuchos::RCP > +integratorForwardSensitivity() +{ + Teuchos::RCP > integrator = + Teuchos::rcp(new Tempus::IntegratorForwardSensitivity()); + return(integrator); +} + +} // namespace Tempus +#endif // Tempus_IntegratorForwardSensitivity_impl_hpp diff --git a/packages/tempus/src/Tempus_IntegratorPseudoTransientForwardSensitivity.cpp b/packages/tempus/src/Tempus_IntegratorPseudoTransientForwardSensitivity.cpp new file mode 100644 index 000000000000..1b5a43e9c817 --- /dev/null +++ b/packages/tempus/src/Tempus_IntegratorPseudoTransientForwardSensitivity.cpp @@ -0,0 +1,37 @@ +// @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_IntegratorPseudoTransientForwardSensitivity.hpp" +#include "Tempus_IntegratorPseudoTransientForwardSensitivity_impl.hpp" + +namespace Tempus { + + TEMPUS_INSTANTIATE_TEMPLATE_CLASS(IntegratorPseudoTransientForwardSensitivity) + + // non-member ctor + template Teuchos::RCP > + integratorPseudoTransientForwardSensitivity( + Teuchos::RCP parameterList, + const Teuchos::RCP >& model); + + // non-member ctor + template Teuchos::RCP > + integratorPseudoTransientForwardSensitivity( + const Teuchos::RCP >& model, + std::string stepperType); + + // non-member ctor + template Teuchos::RCP > + integratorPseudoTransientForwardSensitivity(); + +} // namespace Tempus + +#endif diff --git a/packages/tempus/src/Tempus_IntegratorPseudoTransientForwardSensitivity_decl.hpp b/packages/tempus/src/Tempus_IntegratorPseudoTransientForwardSensitivity_decl.hpp new file mode 100644 index 000000000000..cb82de6114f5 --- /dev/null +++ b/packages/tempus/src/Tempus_IntegratorPseudoTransientForwardSensitivity_decl.hpp @@ -0,0 +1,190 @@ +// @HEADER +// **************************************************************************** +// Tempus: Copyright (2017) Sandia Corporation +// +// Distributed under BSD 3-clause license (See accompanying file Copyright.txt) +// **************************************************************************** +// @HEADER + +#ifndef Tempus_IntegratorPseudoTransientForwardSensitivity_decl_hpp +#define Tempus_IntegratorPseudoTransientForwardSensitivity_decl_hpp + +// Tempus +#include "Tempus_IntegratorBasic.hpp" +#include "Tempus_StaggeredForwardSensitivityModelEvaluator.hpp" + +namespace Tempus { + + +/** \brief Time integrator suitable for pseudotransient forward sensitivity + * analysis */ +/** + * For some problems, time integrators are used to compute steady-state + * solutions (also known as pseudo-transient solvers). When computing + * sensitivities, it is not necessary in these cases to propagate sensitivities + * all the way through the forward time integration. Instead the steady-state + * is first computed as usual, and then the sensitivities are computed using + * a similar pseudo-transient time integration applied to the sensitivity + * equations with the state frozen to the computed steady-state. This + * integrator specializes the transient sensitivity methods implemented by + * Tempus::IntegratorForwardSensitivity to this case. + */ +template +class IntegratorPseudoTransientForwardSensitivity : + virtual public Tempus::Integrator +{ +public: + + /** \brief Constructor with ParameterList and model, and will be fully + * initialized. */ + /*! + * In addition to all of the regular integrator options, the supplied + * parameter list supports the following options contained within a sublist + * "Sensitivities" from the top-level parameter list: + *
    + *
  • "Reuse State Linear Solver" (default: false) Whether to reuse the + * model's W matrix, solver, and preconditioner when solving the + * sensitivity equations. If they can be reused, substantial savings + * in compute time are possible. + *
  • "Force W Update" (default: false) When reusing the solver as above + * whether to force recomputation of W. This can be necessary when + * the solver overwrites it during the solve-phase (e.g., by a + * factorization). + *
  • "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 in the Tempus::CombinedForwardSensitivityModelEvaluator + * documentation. + *
  • "Sensitivity Parameter Index" (default: 0) Model evaluator + * parameter index for which sensitivities will be computed. + *
  • "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. + *
  • "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. + *
  • "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). + *
+ */ + IntegratorPseudoTransientForwardSensitivity( + Teuchos::RCP pList, + const Teuchos::RCP >& model); + + /** \brief Constructor with model and "Stepper Type" and is fully initialized with default settings. */ + IntegratorPseudoTransientForwardSensitivity( + const Teuchos::RCP >& model, + std::string stepperType); + + /// Destructor + /** \brief Constructor that requires a subsequent setParameterList, setStepper, and initialize calls. */ + IntegratorPseudoTransientForwardSensitivity(); + + /// Destructor + virtual ~IntegratorPseudoTransientForwardSensitivity() {} + + /// \name Basic integrator methods + //@{ + + /// Advance the solution to timeMax, and return true if successful. + virtual bool advanceTime(); + /// Advance the solution to timeFinal, and return true if successful. + virtual bool advanceTime(const Scalar timeFinal) override; + /// Get current time + virtual Scalar getTime() const override; + /// Get current index + virtual Scalar getIndex() const override; + /// Get Status + virtual Status getStatus() const override; + /// Get the Stepper + virtual Teuchos::RCP > getStepper() const override; + /// Return a copy of the Tempus ParameterList + virtual Teuchos::RCP getTempusParameterList() override; + virtual void setTempusParameterList(Teuchos::RCP pl) override; + /// Get the SolutionHistory + virtual Teuchos::RCP > getSolutionHistory() const override; + /// Get the TimeStepControl + virtual Teuchos::RCP > getTimeStepControl() const override; + + //@} + + /// Set the initial state from Thyra::VectorBase(s) + virtual void setInitialState( + Scalar t0, + Teuchos::RCP > x0, + Teuchos::RCP > xdot0 = Teuchos::null, + Teuchos::RCP > xdotdot0 = Teuchos::null, + Teuchos::RCP > DxDp0 = Teuchos::null, + Teuchos::RCP > DxdotDp0 = Teuchos::null, + Teuchos::RCP > DxdotdotDp0 = Teuchos::null); + + /// Get current the solution, x + virtual Teuchos::RCP > getX() const; + virtual Teuchos::RCP > getDxDp() const; + /// Get current the time derivative of the solution, xdot + virtual Teuchos::RCP > getXdot() const; + virtual Teuchos::RCP > getDxdotDp() const; + /// Get current the second time derivative of the solution, xdotdot + virtual Teuchos::RCP > getXdotdot() const; + virtual Teuchos::RCP > getDxdotdotDp() const; + + /// \name Overridden from Teuchos::ParameterListAcceptor + //@{ + void setParameterList(const Teuchos::RCP & pl); + Teuchos::RCP getNonconstParameterList() override; + Teuchos::RCP unsetParameterList() override; + + Teuchos::RCP getValidParameters() + const override; + //@} + + /// \name Overridden from Teuchos::Describable + //@{ + std::string description() const override; + void describe(Teuchos::FancyOStream & out, + const Teuchos::EVerbosityLevel verbLevel) const override; + //@} + +protected: + + // Create sensitivity model evaluator from application model + Teuchos::RCP > + createSensitivityModel( + const Teuchos::RCP >& model, + const Teuchos::RCP& inputPL); + + void buildSolutionHistory(); + + Teuchos::RCP > model_; + Teuchos::RCP > sens_model_; + Teuchos::RCP > state_integrator_; + Teuchos::RCP > sens_integrator_; + Teuchos::RCP > solutionHistory_; + bool reuse_solver_; + bool force_W_update_; +}; + +/// Non-member constructor +template +Teuchos::RCP > +integratorPseudoTransientForwardSensitivity( + Teuchos::RCP pList, + const Teuchos::RCP >& model); + +/// Non-member constructor +template +Teuchos::RCP > +integratorPseudoTransientForwardSensitivity( + const Teuchos::RCP >& model, + std::string stepperType); + +/// Non-member constructor +template +Teuchos::RCP > +integratorPseudoTransientForwardSensitivity(); + +} // namespace Tempus + +#endif // Tempus_IntegratorPseudoTransientForwardSensitivity_decl_hpp diff --git a/packages/tempus/src/Tempus_IntegratorPseudoTransientForwardSensitivity_impl.hpp b/packages/tempus/src/Tempus_IntegratorPseudoTransientForwardSensitivity_impl.hpp new file mode 100644 index 000000000000..f4d2a27087ad --- /dev/null +++ b/packages/tempus/src/Tempus_IntegratorPseudoTransientForwardSensitivity_impl.hpp @@ -0,0 +1,581 @@ +// @HEADER +// **************************************************************************** +// Tempus: Copyright (2017) Sandia Corporation +// +// Distributed under BSD 3-clause license (See accompanying file Copyright.txt) +// **************************************************************************** +// @HEADER + +#ifndef Tempus_IntegratorPseudoTransientForwardSensitivity_impl_hpp +#define Tempus_IntegratorPseudoTransientForwardSensitivity_impl_hpp + +#include "Teuchos_VerboseObjectParameterListHelpers.hpp" +#include "Thyra_DefaultMultiVectorProductVector.hpp" +#include "NOX_Thyra.H" + +namespace Tempus { + +template +IntegratorPseudoTransientForwardSensitivity:: +IntegratorPseudoTransientForwardSensitivity( + Teuchos::RCP inputPL, + const Teuchos::RCP >& model) : + reuse_solver_(false) +{ + model_ = model; + sens_model_ = createSensitivityModel(model_, inputPL); + state_integrator_ = integratorBasic(inputPL, model_); + sens_integrator_ = integratorBasic(inputPL, sens_model_); +} + +template +IntegratorPseudoTransientForwardSensitivity:: +IntegratorPseudoTransientForwardSensitivity( + const Teuchos::RCP >& model, + std::string stepperType) : + reuse_solver_(false), + force_W_update_(false) +{ + model_ = model; + sens_model_ = createSensitivityModel(model, Teuchos::null); + state_integrator_ = integratorBasic(model_, stepperType); + sens_integrator_ = integratorBasic(sens_model_, stepperType); +} + +template +IntegratorPseudoTransientForwardSensitivity:: +IntegratorPseudoTransientForwardSensitivity() : + reuse_solver_(false), + force_W_update_(false) +{ + state_integrator_ = integratorBasic(); + sens_integrator_ = integratorBasic(); +} + +template +bool +IntegratorPseudoTransientForwardSensitivity:: +advanceTime() +{ + using Teuchos::RCP; + using Thyra::VectorBase; + + // Run state integrator and get solution + bool state_status = state_integrator_->advanceTime(); + RCP > x = state_integrator_->getX(); + RCP > x_dot = state_integrator_->getXdot(); + RCP > x_dot_dot = state_integrator_->getXdotdot(); + + // Set solution in sensitivity ME + sens_model_->setModelX(x); + sens_model_->setModelXDot(x_dot); + sens_model_->setModelXDotDot(x_dot_dot); + + // Reuse state solver if requested + if (reuse_solver_ && + state_integrator_->getStepper()->getSolver() != Teuchos::null) { + RCP nox_solver = + Teuchos::rcp_dynamic_cast( + state_integrator_->getStepper()->getSolver()); + sens_model_->setLO(nox_solver->get_nonconst_W_op(force_W_update_)); + sens_model_->setPO(nox_solver->get_nonconst_prec_op()); + } + + // Run sensitivity integrator + bool sens_status = sens_integrator_->advanceTime(); + + buildSolutionHistory(); + + return state_status && sens_status; +} + +template +bool +IntegratorPseudoTransientForwardSensitivity:: +advanceTime(const Scalar timeFinal) +{ + using Teuchos::RCP; + using Thyra::VectorBase; + + // Run state integrator and get solution + bool state_status = state_integrator_->advanceTime(timeFinal); + RCP > x = state_integrator_->getX(); + RCP > x_dot = state_integrator_->getXdot(); + RCP > x_dot_dot = state_integrator_->getXdotdot(); + + // Set solution in sensitivity ME + sens_model_->setModelX(x); + sens_model_->setModelXDot(x_dot); + sens_model_->setModelXDotDot(x_dot_dot); + + // Reuse state solver if requested + if (reuse_solver_ && + state_integrator_->getStepper()->getSolver() != Teuchos::null) { + RCP nox_solver = + Teuchos::rcp_dynamic_cast( + state_integrator_->getStepper()->getSolver()); + sens_model_->setLO(nox_solver->get_nonconst_W_op(force_W_update_)); + sens_model_->setPO(nox_solver->get_nonconst_prec_op()); + } + + // Run sensitivity integrator + bool sens_status = sens_integrator_->advanceTime(timeFinal); + + buildSolutionHistory(); + + return state_status && sens_status; +} + +template +Scalar +IntegratorPseudoTransientForwardSensitivity:: +getTime() const +{ + return solutionHistory_->getCurrentTime(); +} + +template +Scalar +IntegratorPseudoTransientForwardSensitivity:: +getIndex() const +{ + return solutionHistory_->getCurrentIndex(); +} + +template +Status +IntegratorPseudoTransientForwardSensitivity:: +getStatus() const +{ + Status state_status = state_integrator_->getStatus(); + Status sens_status = sens_integrator_->getStatus(); + if (state_status == FAILED || sens_status == FAILED) + return FAILED; + if (state_status == WORKING || sens_status == WORKING) + return WORKING; + return PASSED; +} + +template +Teuchos::RCP > +IntegratorPseudoTransientForwardSensitivity:: +getStepper() const +{ + return state_integrator_->getStepper(); +} + +template +Teuchos::RCP +IntegratorPseudoTransientForwardSensitivity:: +getTempusParameterList() +{ + return state_integrator_->getTempusParameterList(); +} + +template +void +IntegratorPseudoTransientForwardSensitivity:: +setTempusParameterList(Teuchos::RCP pl) +{ + state_integrator_->setTempusParameterList(pl); + sens_integrator_->setTempusParameterList(pl); +} + +template +Teuchos::RCP > +IntegratorPseudoTransientForwardSensitivity:: +getSolutionHistory() const +{ + return solutionHistory_; +} + +template +Teuchos::RCP > +IntegratorPseudoTransientForwardSensitivity:: +getTimeStepControl() const +{ + return state_integrator_->getTimeStepControl(); +} + +template +void IntegratorPseudoTransientForwardSensitivity:: +setInitialState(Scalar t0, + Teuchos::RCP > x0, + Teuchos::RCP > xdot0, + Teuchos::RCP > xdotdot0, + Teuchos::RCP > DxDp0, + Teuchos::RCP > DxdotDp0, + Teuchos::RCP > DxdotdotDp0) +{ + using Teuchos::RCP; + using Teuchos::rcp_dynamic_cast; + using Thyra::VectorSpaceBase; + using Thyra::assign; + using Thyra::createMember; + typedef Thyra::DefaultMultiVectorProductVector DMVPV; + + // + // Create and initialize product X, Xdot, Xdotdot + + RCP< const VectorSpaceBase > space = sens_model_->get_x_space(); + RCP X = rcp_dynamic_cast(createMember(space)); + RCP Xdot = rcp_dynamic_cast(createMember(space)); + RCP Xdotdot = rcp_dynamic_cast(createMember(space)); + const Scalar zero = Teuchos::ScalarTraits::zero(); + + // x + if (DxDp0 == Teuchos::null) + assign(X->getNonconstMultiVector().ptr(), zero); + else + assign(X->getNonconstMultiVector().ptr(), *DxDp0); + + // xdot + if (DxdotDp0 == Teuchos::null) + assign(Xdot->getNonconstMultiVector().ptr(), zero); + else + assign(Xdot->getNonconstMultiVector().ptr(), *DxdotDp0); + + // xdotdot + if (DxdotDp0 == Teuchos::null) + assign(Xdotdot->getNonconstMultiVector().ptr(), zero); + else + assign(Xdotdot->getNonconstMultiVector().ptr(), *DxdotdotDp0); + + state_integrator_->setInitialState(t0, x0, xdot0, xdotdot0); + sens_integrator_->setInitialState(t0, X, Xdot, Xdotdot); +} + +template +Teuchos::RCP > +IntegratorPseudoTransientForwardSensitivity:: +getX() const +{ + using Teuchos::RCP; + using Teuchos::rcp_dynamic_cast; + typedef Thyra::DefaultMultiVectorProductVector DMVPV; + + RCP X = + rcp_dynamic_cast(solutionHistory_->getCurrentState()->getX()); + return X->getMultiVector()->col(0); +} + +template +Teuchos::RCP > +IntegratorPseudoTransientForwardSensitivity:: +getDxDp() const +{ + using Teuchos::RCP; + using Teuchos::rcp_dynamic_cast; + typedef Thyra::DefaultMultiVectorProductVector DMVPV; + + RCP X = + rcp_dynamic_cast(solutionHistory_->getCurrentState()->getX()); + const int num_param = X->getMultiVector()->domain()->dim()-1; + const Teuchos::Range1D rng(1,num_param); + return X->getMultiVector()->subView(rng); +} + +template +Teuchos::RCP > +IntegratorPseudoTransientForwardSensitivity:: +getXdot() const +{ + using Teuchos::RCP; + using Teuchos::rcp_dynamic_cast; + typedef Thyra::DefaultMultiVectorProductVector DMVPV; + + RCP Xdot = + rcp_dynamic_cast(solutionHistory_->getCurrentState()->getXDot()); + return Xdot->getMultiVector()->col(0); +} + +template +Teuchos::RCP > +IntegratorPseudoTransientForwardSensitivity:: +getDxdotDp() const +{ + using Teuchos::RCP; + using Teuchos::rcp_dynamic_cast; + typedef Thyra::DefaultMultiVectorProductVector DMVPV; + + RCP Xdot = + rcp_dynamic_cast(solutionHistory_->getCurrentState()->getXDot()); + const int num_param = Xdot->getMultiVector()->domain()->dim()-1; + const Teuchos::Range1D rng(1,num_param); + return Xdot->getMultiVector()->subView(rng); +} + +template +Teuchos::RCP > +IntegratorPseudoTransientForwardSensitivity:: +getXdotdot() const +{ + using Teuchos::RCP; + using Teuchos::rcp_dynamic_cast; + typedef Thyra::DefaultMultiVectorProductVector DMVPV; + + RCP Xdotdot = + rcp_dynamic_cast(solutionHistory_->getCurrentState()->getXDotDot()); + return Xdotdot->getMultiVector()->col(0); +} + +template +Teuchos::RCP > +IntegratorPseudoTransientForwardSensitivity:: +getDxdotdotDp() const +{ + using Teuchos::RCP; + using Teuchos::rcp_dynamic_cast; + typedef Thyra::DefaultMultiVectorProductVector DMVPV; + + RCP Xdotdot = + rcp_dynamic_cast(solutionHistory_->getCurrentState()->getXDotDot()); + const int num_param = Xdotdot->getMultiVector()->domain()->dim()-1; + const Teuchos::Range1D rng(1,num_param); + return Xdotdot->getMultiVector()->subView(rng); +} + +template +std::string +IntegratorPseudoTransientForwardSensitivity:: +description() const +{ + std::string name = "Tempus::IntegratorPseudoTransientForwardSensitivity"; + return(name); +} + +template +void +IntegratorPseudoTransientForwardSensitivity:: +describe( + Teuchos::FancyOStream &out, + const Teuchos::EVerbosityLevel verbLevel) const +{ + out << description() << "::describe" << std::endl; + state_integrator_->describe(out, verbLevel); + sens_integrator_->describe(out, verbLevel); +} + +template +void +IntegratorPseudoTransientForwardSensitivity:: +setParameterList(const Teuchos::RCP & inputPL) +{ + state_integrator_->setParameterList(inputPL); + sens_integrator_->setParameterList(inputPL); + reuse_solver_ = + inputPL->sublist("Sensitivities").get("Reuse State Linear Solver", false); + force_W_update_ = + inputPL->sublist("Sensitivities").get("Force W Update", false); +} + +template +Teuchos::RCP +IntegratorPseudoTransientForwardSensitivity:: +unsetParameterList() +{ + state_integrator_->unsetParameterList(); + return sens_integrator_->unsetParameterList(); +} + +template +Teuchos::RCP +IntegratorPseudoTransientForwardSensitivity:: +getValidParameters() const +{ + Teuchos::RCP pl = + Teuchos::rcp(new Teuchos::ParameterList); + Teuchos::RCP integrator_pl = + state_integrator_->getValidParameters(); + Teuchos::RCP sensitivity_pl = + StaggeredForwardSensitivityModelEvaluator::getValidParameters(); + pl->setParameters(*integrator_pl); + pl->sublist("Sensitivities").setParameters(*sensitivity_pl); + pl->sublist("Sensitivities").set("Reuse State Linear Solver", false); + pl->sublist("Sensitivities").set("Force W Update", false); + + return pl; +} + +template +Teuchos::RCP +IntegratorPseudoTransientForwardSensitivity:: +getNonconstParameterList() +{ + return state_integrator_->getNonconstParameterList(); +} + +template +Teuchos::RCP > +IntegratorPseudoTransientForwardSensitivity:: +createSensitivityModel( + const Teuchos::RCP >& model, + const Teuchos::RCP& inputPL) +{ + using Teuchos::rcp; + + Teuchos::RCP pl = Teuchos::parameterList(); + if (inputPL != Teuchos::null) { + *pl = inputPL->sublist("Sensitivities"); + } + reuse_solver_ = pl->get("Reuse State Linear Solver", false); + force_W_update_ = pl->get("Force W Update", true); + pl->remove("Reuse State Linear Solver"); + pl->remove("Force W Update"); + return rcp(new StaggeredForwardSensitivityModelEvaluator(model, pl)); +} + +template +void +IntegratorPseudoTransientForwardSensitivity:: +buildSolutionHistory() +{ + 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::DefaultMultiVectorProductVector DMVPV; + + // Create combined solution histories, first for the states with zero + // sensitivities and then for the sensitivities with frozen states + RCP shPL = + Teuchos::sublist(state_integrator_->getIntegratorParameterList(), + "Solution History", true); + solutionHistory_ = rcp(new SolutionHistory(shPL)); + + const int num_param = + rcp_dynamic_cast(sens_integrator_->getX())->getMultiVector()->domain()->dim(); + RCP > x_space = model_->get_x_space(); + RCP > prod_space = + Thyra::multiVectorProductVectorSpace(x_space, num_param+1); + const Teuchos::Range1D rng(1,num_param); + const Scalar zero = Teuchos::ScalarTraits::zero(); + + RCP > state_solution_history = + state_integrator_->getSolutionHistory(); + int num_states = state_solution_history->getNumStates(); + for (int i=0; i > state = (*state_solution_history)[i]; + + // X + RCP< MultiVectorBase > x_mv = + createMembers(x_space, num_param+1); + assign(x_mv->col(0).ptr(), *(state->getX())); + assign(x_mv->subView(rng).ptr(), zero); + RCP x = multiVectorProductVector(prod_space, x_mv); + + // X-Dot + RCP< MultiVectorBase > x_dot_mv = + createMembers(x_space, num_param+1); + assign(x_dot_mv->col(0).ptr(), *(state->getXDot())); + assign(x_dot_mv->subView(rng).ptr(), zero); + RCP x_dot = multiVectorProductVector(prod_space, x_dot_mv); + + // X-Dot-Dot + RCP< MultiVectorBase > x_dot_dot_mv = + createMembers(x_space, num_param+1); + assign(x_dot_dot_mv->col(0).ptr(), *(state->getXDotDot())); + assign(x_dot_dot_mv->subView(rng).ptr(), zero); + RCP x_dot_dot = multiVectorProductVector(prod_space, x_dot_dot_mv); + + RCP > prod_state = + rcp(new SolutionState(state->getMetaData()->clone(), + x, x_dot, x_dot_dot, + state->getStepperState()->clone())); + solutionHistory_->addState(prod_state); + } + + RCP > frozen_x = + state_solution_history->getCurrentState()->getX(); + RCP > frozen_x_dot = + state_solution_history->getCurrentState()->getXDot(); + RCP > frozen_x_dot_dot = + state_solution_history->getCurrentState()->getXDotDot(); + RCP > sens_solution_history = + sens_integrator_->getSolutionHistory(); + num_states = sens_solution_history->getNumStates(); + for (int i=0; i > state = (*sens_solution_history)[i]; + + // X + RCP< MultiVectorBase > x_mv = + createMembers(x_space, num_param+1); + RCP > dxdp = + rcp_dynamic_cast(state->getX())->getMultiVector(); + assign(x_mv->col(0).ptr(), *(frozen_x)); + assign(x_mv->subView(rng).ptr(), *dxdp); + RCP x = multiVectorProductVector(prod_space, x_mv); + + // X-Dot + RCP< MultiVectorBase > x_dot_mv = + createMembers(x_space, num_param+1); + RCP > dxdotdp = + rcp_dynamic_cast(state->getXDot())->getMultiVector(); + assign(x_dot_mv->col(0).ptr(), *(frozen_x_dot)); + assign(x_dot_mv->subView(rng).ptr(), *dxdotdp); + RCP x_dot = multiVectorProductVector(prod_space, x_dot_mv); + + // X-Dot-Dot + RCP< MultiVectorBase > x_dot_dot_mv = + createMembers(x_space, num_param+1); + RCP > dxdotdotdp = + rcp_dynamic_cast(state->getXDotDot())->getMultiVector(); + assign(x_dot_dot_mv->col(0).ptr(), *(frozen_x_dot_dot)); + assign(x_dot_dot_mv->subView(rng).ptr(), *dxdotdotdp); + RCP x_dot_dot = multiVectorProductVector(prod_space, x_dot_dot_mv); + + RCP > prod_state = + rcp(new SolutionState(state->getMetaData()->clone(), + x, x_dot, x_dot_dot, + state->getStepperState()->clone())); + solutionHistory_->addState(prod_state); + } + + // Set current state to the last added state + RCP > last_state = + (*solutionHistory_)[solutionHistory_->getNumStates()-1]; + solutionHistory_->setCurrentState(last_state); +} + +/// Non-member constructor +template +Teuchos::RCP > +integratorPseudoTransientForwardSensitivity( + Teuchos::RCP pList, + const Teuchos::RCP >& model) +{ + Teuchos::RCP > integrator = + Teuchos::rcp(new Tempus::IntegratorPseudoTransientForwardSensitivity(pList, model)); + return(integrator); +} + +/// Non-member constructor +template +Teuchos::RCP > +integratorPseudoTransientForwardSensitivity( + const Teuchos::RCP >& model, + std::string stepperType) +{ + Teuchos::RCP > integrator = + Teuchos::rcp(new Tempus::IntegratorPseudoTransientForwardSensitivity(model, stepperType)); + return(integrator); +} + +/// Non-member constructor +template +Teuchos::RCP > +integratorPseudoTransientForwardSensitivity() +{ + Teuchos::RCP > integrator = + Teuchos::rcp(new Tempus::IntegratorPseudoTransientForwardSensitivity()); + return(integrator); +} + +} // namespace Tempus +#endif // Tempus_IntegratorPseudoTransientForwardSensitivity_impl_hpp diff --git a/packages/tempus/src/Tempus_SolutionHistory_decl.hpp b/packages/tempus/src/Tempus_SolutionHistory_decl.hpp index bb3d0d1b1368..44d4fb8efad7 100644 --- a/packages/tempus/src/Tempus_SolutionHistory_decl.hpp +++ b/packages/tempus/src/Tempus_SolutionHistory_decl.hpp @@ -120,6 +120,10 @@ class SolutionHistory Teuchos::RCP > getCurrentState() const { return currentState_; } + /// Set the current state, i.e., the last accepted state + void setCurrentState(const Teuchos::RCP >& state) + { currentState_ = state; } + /// Return the working state Teuchos::RCP > getWorkingState() const { return workingState_; } diff --git a/packages/tempus/src/Tempus_SolutionStateMetaData_decl.hpp b/packages/tempus/src/Tempus_SolutionStateMetaData_decl.hpp index 6c064c42d902..e6f16312c32c 100644 --- a/packages/tempus/src/Tempus_SolutionStateMetaData_decl.hpp +++ b/packages/tempus/src/Tempus_SolutionStateMetaData_decl.hpp @@ -51,7 +51,7 @@ class SolutionStateMetaData : SolutionStateMetaData(const SolutionStateMetaData& ssmd); /// Clone constructor - Teuchos::RCP > clone(); + Teuchos::RCP > clone() const; /// This is a deep copy void copy(Teuchos::RCP > ssmd); diff --git a/packages/tempus/src/Tempus_SolutionStateMetaData_impl.hpp b/packages/tempus/src/Tempus_SolutionStateMetaData_impl.hpp index 8c4cceed58b9..0724efe5f6fc 100644 --- a/packages/tempus/src/Tempus_SolutionStateMetaData_impl.hpp +++ b/packages/tempus/src/Tempus_SolutionStateMetaData_impl.hpp @@ -84,7 +84,7 @@ SolutionStateMetaData::SolutionStateMetaData(const SolutionStateMetaData template -Teuchos::RCP > SolutionStateMetaData::clone() +Teuchos::RCP > SolutionStateMetaData::clone() const { Teuchos::RCP > md = rcp(new SolutionStateMetaData ( diff --git a/packages/tempus/src/Tempus_SolutionState_decl.hpp b/packages/tempus/src/Tempus_SolutionState_decl.hpp index 90b24dab8f13..3cdcfa69dec2 100644 --- a/packages/tempus/src/Tempus_SolutionState_decl.hpp +++ b/packages/tempus/src/Tempus_SolutionState_decl.hpp @@ -122,12 +122,16 @@ class SolutionState : /// Return the Stepper status virtual Status getStepperStatus() const - {return stepperState_->stepperStatus_;}; + {return stepperState_->stepperStatus_;} /// Return Meta Data virtual Teuchos::RCP > getMetaData() { return metaData_; } + /// Return Meta Data + virtual Teuchos::RCP > getMetaData() const + { return metaData_; } + /// Get the current solution, x. virtual Teuchos::RCP > getX() {return x_;} @@ -154,6 +158,10 @@ class SolutionState : virtual Teuchos::RCP > getStepperState() { return stepperState_; } + /// Get the StepperState + virtual Teuchos::RCP > getStepperState() const + { return stepperState_; } + //@} diff --git a/packages/tempus/src/Tempus_StaggeredForwardSensitivityModelEvaluator.cpp b/packages/tempus/src/Tempus_StaggeredForwardSensitivityModelEvaluator.cpp new file mode 100644 index 000000000000..34bd94ce7314 --- /dev/null +++ b/packages/tempus/src/Tempus_StaggeredForwardSensitivityModelEvaluator.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_StaggeredForwardSensitivityModelEvaluator_decl.hpp" +#include "Tempus_StaggeredForwardSensitivityModelEvaluator_impl.hpp" + +namespace Tempus { + TEMPUS_INSTANTIATE_TEMPLATE_CLASS(StaggeredForwardSensitivityModelEvaluator) +} + +#endif diff --git a/packages/tempus/src/Tempus_StaggeredForwardSensitivityModelEvaluator_decl.hpp b/packages/tempus/src/Tempus_StaggeredForwardSensitivityModelEvaluator_decl.hpp new file mode 100644 index 000000000000..b098b24a7a49 --- /dev/null +++ b/packages/tempus/src/Tempus_StaggeredForwardSensitivityModelEvaluator_decl.hpp @@ -0,0 +1,168 @@ +// @HEADER +// **************************************************************************** +// Tempus: Copyright (2017) Sandia Corporation +// +// Distributed under BSD 3-clause license (See accompanying file Copyright.txt) +// **************************************************************************** +// @HEADER + +#ifndef Tempus_StaggeredForwardSensitivityModelEvaluator_decl_hpp +#define Tempus_StaggeredForwardSensitivityModelEvaluator_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 the forward sensitivity equations: + * F(X) = (df/dx)(x,p) * X + df/dp(x,p) = 0 + * where 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). + * + * This model evaluator differes from CombinedForwardSensitivityModelEvaluator + * in that it doesn't include the state equations in the residual, just the + * sensitivity equations. Therefore it provides methods to set the state + * solution vector (x) and time derivatives (x_dot, x_dot_dot) for use in + * evaluating the senstivity residual. It also provides methods for setting + * the linear operator (W a.k.a. alpha*df/dx + beta*df/dx_dot) and its + * preconditioner in cases where they can be reused from the state model + * evaluations. + */ +template +class StaggeredForwardSensitivityModelEvaluator : + public Thyra::StateFuncModelEvaluatorBase { +public: + typedef Thyra::VectorBase Vector; + typedef Thyra::MultiVectorBase MultiVector; + + //! Constructor + /*! + * The optionally supplied parameter list supports the following options: + *
    + *
  • "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(). + *
  • "Sensitivity Parameter Index" (default: 0) Model evaluator + * parameter index for which sensitivities will be computed. + *
  • "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. + *
  • "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. + *
  • "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). + *
+ */ + StaggeredForwardSensitivityModelEvaluator( + const Teuchos::RCP > & model, + const Teuchos::RCP& pList = Teuchos::null, + const Teuchos::RCP& dxdp_init = Teuchos::null, + const Teuchos::RCP& dx_dotdp_init = Teuchos::null, + const Teuchos::RCP& dx_dotdot_dp_init = Teuchos::null); + + //! Get the underlying model 'f' + Teuchos::RCP > getModel() const + { return model_; } + + //! Set the x vector corresponding to the underlying model + void setModelX(const Teuchos::RCP >& x) + { x_ = x; } + + //! Set the x_dot vector corresponding to the underlying model + void setModelXDot(const Teuchos::RCP >& x_dot) + { x_dot_ = x_dot; } + + //! Set the x_dot_dot vector corresponding to the underlying model + void setModelXDotDot(const Teuchos::RCP >& x_dot_dot) + { x_dot_dot_ = x_dot_dot; } + + //! Set the LO (W) of the underlying model if you want to reuse it + void setLO(const Teuchos::RCP >& lo) + { lo_ = lo; } + + //! Set the preconditioner of the underlying model if you want to reuse it + void setPO(const Teuchos::RCP >& po) + { po_ = po; } + + /** \name Public functions overridden from ModelEvaulator. */ + //@{ + + Teuchos::RCP > get_p_space(int p) const; + + Teuchos::RCP > get_p_names(int p) const; + + Teuchos::RCP > get_x_space() const; + + Teuchos::RCP > get_f_space() const; + + Teuchos::RCP > create_W_op() const; + + Teuchos::RCP > + get_W_factory() const; + + Thyra::ModelEvaluatorBase::InArgs createInArgs() const; + + Thyra::ModelEvaluatorBase::InArgs getNominalValues() const; + + //@} + + static Teuchos::RCP getValidParameters(); + +private: + + Thyra::ModelEvaluatorBase::OutArgs createOutArgsImpl() const; + + void evalModelImpl( + const Thyra::ModelEvaluatorBase::InArgs &inArgs, + const Thyra::ModelEvaluatorBase::OutArgs &outArgs) const; + + + Thyra::ModelEvaluatorBase::InArgs prototypeInArgs_; + Thyra::ModelEvaluatorBase::OutArgs prototypeOutArgs_; + + Teuchos::RCP > model_; + Teuchos::RCP dxdp_init_; + Teuchos::RCP dx_dotdp_init_; + Teuchos::RCP 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 > dxdp_space_; + Teuchos::RCP > dfdp_space_; + Teuchos::RCP > x_; + Teuchos::RCP > x_dot_; + Teuchos::RCP > x_dot_dot_; + + Teuchos::RCP > lo_; + Teuchos::RCP > po_; + + mutable Teuchos::RCP > my_dfdx_; + mutable Teuchos::RCP > my_dfdxdot_; + mutable Teuchos::RCP > my_dfdxdotdot_; +}; + +} // namespace Tempus + +#endif diff --git a/packages/tempus/src/Tempus_StaggeredForwardSensitivityModelEvaluator_impl.hpp b/packages/tempus/src/Tempus_StaggeredForwardSensitivityModelEvaluator_impl.hpp new file mode 100644 index 000000000000..495468dc707e --- /dev/null +++ b/packages/tempus/src/Tempus_StaggeredForwardSensitivityModelEvaluator_impl.hpp @@ -0,0 +1,381 @@ +// @HEADER +// **************************************************************************** +// Tempus: Copyright (2017) Sandia Corporation +// +// Distributed under BSD 3-clause license (See accompanying file Copyright.txt) +// **************************************************************************** +// @HEADER + +#ifndef Tempus_StaggeredForwardSensitivityModelEvaluator_impl_hpp +#define Tempus_StaggeredForwardSensitivityModelEvaluator_impl_hpp + +#include "Thyra_DefaultMultiVectorProductVector.hpp" +#include "Thyra_MultiVectorLinearOp.hpp" +#include "Thyra_MultiVectorLinearOpWithSolveFactory.hpp" +#include "Thyra_ReuseLinearOpWithSolveFactory.hpp" + +namespace Tempus { + +template +StaggeredForwardSensitivityModelEvaluator:: +StaggeredForwardSensitivityModelEvaluator( + const Teuchos::RCP > & model, + const Teuchos::RCP& pList, + const Teuchos::RCP& dxdp_init, + const Teuchos::RCP& dx_dotdp_init, + const Teuchos::RCP& dx_dotdotdp_init) : + model_(model), + dxdp_init_(dxdp_init), + dx_dotdp_init_(dx_dotdp_init), + dx_dotdotdp_init_(dx_dotdotdp_init), + p_index_(0), + x_tangent_index_(1), + xdot_tangent_index_(2), + xdotdot_tangent_index_(3), + use_dfdp_as_tangent_(false) +{ + typedef Thyra::ModelEvaluatorBase MEB; + + // Set parameters + Teuchos::RCP pl = + Teuchos::rcp(new Teuchos::ParameterList); + if (pList != Teuchos::null) + *pl = *pList; + pl->validateParametersAndSetDefaults(*this->getValidParameters()); + use_dfdp_as_tangent_ = pl->get("Use DfDp as Tangent"); + p_index_ = pl->get("Sensitivity Parameter Index"); + x_tangent_index_ = pl->get("Sensitivity X Tangent Index"); + xdot_tangent_index_ = pl->get("Sensitivity X-Dot Tangent Index"); + xdotdot_tangent_index_ = pl->get("Sensitivity X-Dot-Dot Tangent Index"); + + num_param_ = model_->get_p_space(p_index_)->dim(); + dxdp_space_ = + Thyra::multiVectorProductVectorSpace(model_->get_x_space(), num_param_); + dfdp_space_ = + Thyra::multiVectorProductVectorSpace(model_->get_f_space(), num_param_); + + MEB::InArgs me_inArgs = model_->createInArgs(); + MEB::InArgsSetup inArgs; + inArgs.setModelEvalDescription(this->description()); + inArgs.setSupports(MEB::IN_ARG_x); + if (me_inArgs.supports(MEB::IN_ARG_x_dot)) + inArgs.setSupports(MEB::IN_ARG_x_dot); + if (me_inArgs.supports(MEB::IN_ARG_t)) + inArgs.setSupports(MEB::IN_ARG_t); + if (me_inArgs.supports(MEB::IN_ARG_alpha)) + inArgs.setSupports(MEB::IN_ARG_alpha); + if (me_inArgs.supports(MEB::IN_ARG_beta)) + inArgs.setSupports(MEB::IN_ARG_beta); + if (me_inArgs.supports(MEB::IN_ARG_W_x_dot_dot_coeff)) + inArgs.setSupports(MEB::IN_ARG_W_x_dot_dot_coeff); + + // Support additional parameters for x and xdot + inArgs.set_Np(me_inArgs.Np()); + prototypeInArgs_ = inArgs; + + MEB::OutArgs me_outArgs = model_->createOutArgs(); + MEB::OutArgsSetup outArgs; + outArgs.setModelEvalDescription(this->description()); + outArgs.set_Np_Ng(me_inArgs.Np(),0); + outArgs.setSupports(MEB::OUT_ARG_f); + outArgs.setSupports(MEB::OUT_ARG_W_op); + prototypeOutArgs_ = outArgs; + + TEUCHOS_ASSERT(me_outArgs.supports(MEB::OUT_ARG_DfDp, p_index_).supports(MEB::DERIV_MV_JACOBIAN_FORM)); + if (!use_dfdp_as_tangent_) + TEUCHOS_ASSERT(me_outArgs.supports(MEB::OUT_ARG_W_op)); +} + +template +Teuchos::RCP > +StaggeredForwardSensitivityModelEvaluator:: +get_p_space(int p) const +{ + TEUCHOS_ASSERT(p < model_->Np()); + return model_->get_p_space(p); +} + +template +Teuchos::RCP > +StaggeredForwardSensitivityModelEvaluator:: +get_p_names(int p) const +{ + TEUCHOS_ASSERT(p < model_->Np()); + return model_->get_p_names(p); +} + +template +Teuchos::RCP > +StaggeredForwardSensitivityModelEvaluator:: +get_x_space() const +{ + return dxdp_space_; +} + +template +Teuchos::RCP > +StaggeredForwardSensitivityModelEvaluator:: +get_f_space() const +{ + return dfdp_space_; +} + +template +Teuchos::RCP > +StaggeredForwardSensitivityModelEvaluator:: +create_W_op() const +{ + Teuchos::RCP > op; + if (lo_ != Teuchos::null) + op = lo_; + else + op = model_->create_W_op(); + return Thyra::nonconstMultiVectorLinearOp(op, num_param_); +} + +template +Teuchos::RCP > +StaggeredForwardSensitivityModelEvaluator:: +get_W_factory() const +{ + typedef Thyra::DefaultMultiVectorProductVectorSpace DMVPVS; + using Teuchos::RCP; + using Teuchos::rcp_dynamic_cast; + + Teuchos::RCP > factory; + Teuchos::RCP > model_factory + = model_->get_W_factory(); + if (po_ != Teuchos::null) { + factory = + Thyra::reuseLinearOpWithSolveFactory(model_factory, po_); + } + else + factory = model_factory; + RCP > mv_op = this->create_W_op(); + RCP mv_domain = + rcp_dynamic_cast(mv_op->domain()); + RCP mv_range = + rcp_dynamic_cast(mv_op->range()); + return Thyra::multiVectorLinearOpWithSolveFactory(factory, + mv_range, + mv_domain); +} + +template +Thyra::ModelEvaluatorBase::InArgs +StaggeredForwardSensitivityModelEvaluator:: +createInArgs() const +{ + return prototypeInArgs_; +} + +template +Thyra::ModelEvaluatorBase::InArgs +StaggeredForwardSensitivityModelEvaluator:: +getNominalValues() const +{ + typedef Thyra::ModelEvaluatorBase MEB; + typedef Thyra::DefaultMultiVectorProductVector DMVPV; + using Teuchos::RCP; + using Teuchos::rcp_dynamic_cast; + + MEB::InArgs me_nominal = model_->getNominalValues(); + MEB::InArgs nominal = this->createInArgs(); + + const Scalar zero = Teuchos::ScalarTraits::zero(); + + // Set initial x. If dxdp_init == null, set initial dx/dp = 0 + RCP< Thyra::VectorBase > x = Thyra::createMember(*dxdp_space_); + RCP dxdp = rcp_dynamic_cast(x); + if (dxdp_init_ == Teuchos::null) + Thyra::assign(dxdp->getNonconstMultiVector().ptr(), zero); + else + Thyra::assign(dxdp->getNonconstMultiVector().ptr(), *dxdp_init_); + nominal.set_x(x); + + // Set initial xdot. If dx_dotdp_init == null, set initial dxdot/dp = 0 + if (me_nominal.supports(MEB::IN_ARG_x_dot)) { + RCP< Thyra::VectorBase > xdot = Thyra::createMember(*dxdp_space_); + RCP dxdotdp = rcp_dynamic_cast(xdot); + if (dx_dotdp_init_ == Teuchos::null) + Thyra::assign(dxdotdp->getNonconstMultiVector().ptr(), zero); + else + Thyra::assign(dxdotdp->getNonconstMultiVector().ptr(), *dx_dotdp_init_); + nominal.set_x_dot(xdot); + } + + // Set initial xdotdot. If dx_dotdotdp_init == null, set initial dxdotdot/dp = 0 + if (me_nominal.supports(MEB::IN_ARG_x_dot_dot)) { + RCP< Thyra::VectorBase > xdotdot = + Thyra::createMember(*dxdp_space_); + RCP dxdotdotdp = rcp_dynamic_cast(xdotdot); + if (dx_dotdotdp_init_ == Teuchos::null) + Thyra::assign(dxdotdotdp->getNonconstMultiVector().ptr(), zero); + else + Thyra::assign(dxdotdotdp->getNonconstMultiVector().ptr(), + *dx_dotdotdp_init_); + nominal.set_x_dot_dot(xdotdot); + } + + const int np = model_->Np(); + for (int i=0; i +Thyra::ModelEvaluatorBase::OutArgs +StaggeredForwardSensitivityModelEvaluator:: +createOutArgsImpl() const +{ + return prototypeOutArgs_; +} + +template +void +StaggeredForwardSensitivityModelEvaluator:: +evalModelImpl(const Thyra::ModelEvaluatorBase::InArgs &inArgs, + const Thyra::ModelEvaluatorBase::OutArgs &outArgs) const +{ + typedef Thyra::ModelEvaluatorBase MEB; + typedef Thyra::DefaultMultiVectorProductVector DMVPV; + using Teuchos::RCP; + using Teuchos::rcp_dynamic_cast; + + // setup input arguments for model + RCP< const Thyra::MultiVectorBase > dxdp, dxdotdp, dxdotdotdp; + MEB::InArgs me_inArgs = model_->getNominalValues(); + dxdp = rcp_dynamic_cast(inArgs.get_x())->getMultiVector(); + TEUCHOS_ASSERT(x_ != Teuchos::null); + me_inArgs.set_x(x_); + if (use_dfdp_as_tangent_) + me_inArgs.set_p(x_tangent_index_, inArgs.get_x()); + if (me_inArgs.supports(MEB::IN_ARG_x_dot)) { + if (inArgs.get_x_dot() != Teuchos::null) { + dxdotdp = + rcp_dynamic_cast(inArgs.get_x_dot())->getMultiVector(); + TEUCHOS_ASSERT(x_dot_ != Teuchos::null); + me_inArgs.set_x_dot(x_dot_); + if (use_dfdp_as_tangent_) + me_inArgs.set_p(xdot_tangent_index_, inArgs.get_x_dot()); + } + else // clear out xdot if it was set in nominalValues + me_inArgs.set_x_dot(Teuchos::null); + } + if (me_inArgs.supports(MEB::IN_ARG_x_dot_dot)) { + if (inArgs.get_x_dot_dot() != Teuchos::null) { + dxdotdotdp = + rcp_dynamic_cast(inArgs.get_x_dot_dot())->getMultiVector(); + TEUCHOS_ASSERT(x_dot_dot_ != Teuchos::null); + me_inArgs.set_x_dot_dot(x_dot_dot_); + if (use_dfdp_as_tangent_) + me_inArgs.set_p(xdotdot_tangent_index_, inArgs.get_x_dot_dot()); + } + else // clear out xdotdot if it was set in nominalValues + me_inArgs.set_x_dot_dot(Teuchos::null); + } + if (me_inArgs.supports(MEB::IN_ARG_t)) + me_inArgs.set_t(inArgs.get_t()); + if (me_inArgs.supports(MEB::IN_ARG_alpha)) + me_inArgs.set_alpha(inArgs.get_alpha()); + if (me_inArgs.supports(MEB::IN_ARG_beta)) + me_inArgs.set_beta(inArgs.get_beta()); + if (me_inArgs.supports(MEB::IN_ARG_W_x_dot_dot_coeff)) + me_inArgs.set_W_x_dot_dot_coeff(inArgs.get_W_x_dot_dot_coeff()); + + // Set parameters -- be careful to only set ones that were set in our + // inArgs to not null out any specified through nominalValues or + // dx/dp above + const int np = me_inArgs.Np(); + for (int i=0; i > dfdp; + MEB::OutArgs me_outArgs = model_->createOutArgs(); + if (outArgs.get_f() != Teuchos::null) { + dfdp = rcp_dynamic_cast(outArgs.get_f())->getNonconstMultiVector(); + me_outArgs.set_DfDp(p_index_, dfdp); + } + if (lo_ == Teuchos::null && outArgs.get_W_op() != Teuchos::null) { + RCP > op = outArgs.get_W_op(); + RCP > mv_op = + rcp_dynamic_cast >(op); + me_outArgs.set_W_op(mv_op->getNonconstLinearOp()); + } + + // build residual and jacobian + model_->evalModel(me_inArgs, me_outArgs); + + // Compute (df/dx) * (dx/dp) + (df/dxdot) * (dxdot/dp) + (df/dxdotdot) * (dxdotdot/dp) + (df/dp) + // if the underlying ME doesn't already do this. This requires computing + // df/dx, df/dxdot, df/dxdotdot as separate operators + if (!use_dfdp_as_tangent_) { + if (dxdp != Teuchos::null && dfdp != Teuchos::null) { + if (my_dfdx_ == Teuchos::null) + my_dfdx_ = model_->create_W_op(); + MEB::OutArgs meo = model_->createOutArgs(); + meo.set_W_op(my_dfdx_); + if (me_inArgs.supports(MEB::IN_ARG_alpha)) + me_inArgs.set_alpha(0.0); + if (me_inArgs.supports(MEB::IN_ARG_beta)) + me_inArgs.set_beta(1.0); + if (me_inArgs.supports(MEB::IN_ARG_W_x_dot_dot_coeff)) + me_inArgs.set_W_x_dot_dot_coeff(0.0); + model_->evalModel(me_inArgs, meo); + my_dfdx_->apply(Thyra::NOTRANS, *dxdp, dfdp.ptr(), Scalar(1.0), Scalar(1.0)); + } + if (dxdotdp != Teuchos::null && dfdp != Teuchos::null) { + if (my_dfdxdot_ == Teuchos::null) + my_dfdxdot_ = model_->create_W_op(); + MEB::OutArgs meo = model_->createOutArgs(); + meo.set_W_op(my_dfdxdot_); + if (me_inArgs.supports(MEB::IN_ARG_alpha)) + me_inArgs.set_alpha(1.0); + if (me_inArgs.supports(MEB::IN_ARG_beta)) + me_inArgs.set_beta(0.0); + if (me_inArgs.supports(MEB::IN_ARG_W_x_dot_dot_coeff)) + me_inArgs.set_W_x_dot_dot_coeff(0.0); + model_->evalModel(me_inArgs, meo); + my_dfdxdot_->apply(Thyra::NOTRANS, *dxdotdp, dfdp.ptr(), Scalar(1.0), Scalar(1.0)); + } + if (dxdotdotdp != Teuchos::null && dfdp != Teuchos::null) { + if (my_dfdxdotdot_ == Teuchos::null) + my_dfdxdotdot_ = model_->create_W_op(); + MEB::OutArgs meo = model_->createOutArgs(); + meo.set_W_op(my_dfdxdotdot_); + if (me_inArgs.supports(MEB::IN_ARG_alpha)) + me_inArgs.set_alpha(0.0); + if (me_inArgs.supports(MEB::IN_ARG_beta)) + me_inArgs.set_beta(0.0); + if (me_inArgs.supports(MEB::IN_ARG_W_x_dot_dot_coeff)) + me_inArgs.set_W_x_dot_dot_coeff(1.0); + model_->evalModel(me_inArgs, meo); + my_dfdxdotdot_->apply(Thyra::NOTRANS, *dxdotdotdp, dfdp.ptr(), Scalar(1.0), Scalar(1.0)); + } + } +} + +template +Teuchos::RCP +StaggeredForwardSensitivityModelEvaluator:: +getValidParameters() +{ + Teuchos::RCP pl = Teuchos::parameterList(); + pl->set("Use DfDp as Tangent", false); + pl->set("Sensitivity Parameter Index", 0); + pl->set("Sensitivity X Tangent Index", 1); + pl->set("Sensitivity X-Dot Tangent Index", 2); + pl->set("Sensitivity X-Dot-Dot Tangent Index", 3); + return pl; +} + +} // namespace Tempus + +#endif diff --git a/packages/tempus/src/Tempus_Stepper.hpp b/packages/tempus/src/Tempus_Stepper.hpp index efa6df4c40a5..6c606464ce99 100644 --- a/packages/tempus/src/Tempus_Stepper.hpp +++ b/packages/tempus/src/Tempus_Stepper.hpp @@ -84,6 +84,8 @@ class Stepper /// Set solver. virtual void setSolver( Teuchos::RCP > solver) = 0; + /// Get solver + virtual Teuchos::RCP > getSolver() const = 0; /// Initialize during construction and after changing input parameters. virtual void initialize() = 0; diff --git a/packages/tempus/src/Tempus_StepperBDF2_decl.hpp b/packages/tempus/src/Tempus_StepperBDF2_decl.hpp index 82e515cf5958..22883148653e 100644 --- a/packages/tempus/src/Tempus_StepperBDF2_decl.hpp +++ b/packages/tempus/src/Tempus_StepperBDF2_decl.hpp @@ -60,6 +60,8 @@ class StepperBDF2 : virtual public Tempus::StepperImplicit Teuchos::RCP solverPL=Teuchos::null); virtual void setSolver( Teuchos::RCP > solver); + virtual Teuchos::RCP > getSolver() const + { return solver_; } virtual void setObserver( Teuchos::RCP > obs = Teuchos::null); diff --git a/packages/tempus/src/Tempus_StepperBackwardEuler_decl.hpp b/packages/tempus/src/Tempus_StepperBackwardEuler_decl.hpp index 87889d1543a9..3c81b146fc5f 100644 --- a/packages/tempus/src/Tempus_StepperBackwardEuler_decl.hpp +++ b/packages/tempus/src/Tempus_StepperBackwardEuler_decl.hpp @@ -52,6 +52,8 @@ class StepperBackwardEuler : virtual public Tempus::StepperImplicit Teuchos::RCP solverPL=Teuchos::null); virtual void setSolver( Teuchos::RCP > solver); + virtual Teuchos::RCP > getSolver() const + { return solver_; } virtual void setObserver( Teuchos::RCP > obs = Teuchos::null); diff --git a/packages/tempus/src/Tempus_StepperDIRK_decl.hpp b/packages/tempus/src/Tempus_StepperDIRK_decl.hpp index 4439ff3197c5..e90893c10183 100644 --- a/packages/tempus/src/Tempus_StepperDIRK_decl.hpp +++ b/packages/tempus/src/Tempus_StepperDIRK_decl.hpp @@ -105,6 +105,8 @@ class StepperDIRK : virtual public Tempus::StepperImplicit Teuchos::RCP solverPL = Teuchos::null); virtual void setSolver( Teuchos::RCP > solver); + virtual Teuchos::RCP > getSolver() const + { return solver_; } virtual void setObserver( Teuchos::RCP > obs = Teuchos::null); diff --git a/packages/tempus/src/Tempus_StepperExplicitRK_decl.hpp b/packages/tempus/src/Tempus_StepperExplicitRK_decl.hpp index 72be57b02e9b..eefe0b9516ce 100644 --- a/packages/tempus/src/Tempus_StepperExplicitRK_decl.hpp +++ b/packages/tempus/src/Tempus_StepperExplicitRK_decl.hpp @@ -85,6 +85,8 @@ class StepperExplicitRK : virtual public Tempus::Stepper Teuchos::RCP solverPL=Teuchos::null); virtual void setSolver( Teuchos::RCP > solver); + virtual Teuchos::RCP > getSolver() const + { return Teuchos::null; } virtual void setObserver( Teuchos::RCP > obs = Teuchos::null); diff --git a/packages/tempus/src/Tempus_StepperForwardEuler_decl.hpp b/packages/tempus/src/Tempus_StepperForwardEuler_decl.hpp index 4e9fae4536a9..104b90921294 100644 --- a/packages/tempus/src/Tempus_StepperForwardEuler_decl.hpp +++ b/packages/tempus/src/Tempus_StepperForwardEuler_decl.hpp @@ -62,6 +62,8 @@ class StepperForwardEuler : virtual public Tempus::Stepper Teuchos::RCP solverPL=Teuchos::null); virtual void setSolver( Teuchos::RCP > solver); + virtual Teuchos::RCP > getSolver() const + { return Teuchos::null; } virtual void setObserver( Teuchos::RCP > obs = Teuchos::null); diff --git a/packages/tempus/src/Tempus_StepperHHTAlpha_decl.hpp b/packages/tempus/src/Tempus_StepperHHTAlpha_decl.hpp index a021f0f85395..2286714c93d6 100644 --- a/packages/tempus/src/Tempus_StepperHHTAlpha_decl.hpp +++ b/packages/tempus/src/Tempus_StepperHHTAlpha_decl.hpp @@ -62,6 +62,8 @@ class StepperHHTAlpha : virtual public Tempus::StepperImplicit Teuchos::RCP solverPL=Teuchos::null); virtual void setSolver( Teuchos::RCP > solver); + virtual Teuchos::RCP > getSolver() const + { return solver_; } /// Initialize during construction and after changing input parameters. diff --git a/packages/tempus/src/Tempus_StepperIMEX_RK_Partition_decl.hpp b/packages/tempus/src/Tempus_StepperIMEX_RK_Partition_decl.hpp index 8cbab8ce2453..ed5e8f8651d6 100644 --- a/packages/tempus/src/Tempus_StepperIMEX_RK_Partition_decl.hpp +++ b/packages/tempus/src/Tempus_StepperIMEX_RK_Partition_decl.hpp @@ -282,6 +282,8 @@ class StepperIMEX_RK_Partition : virtual public Tempus::StepperImplicit /// Set solver. virtual void setSolver( Teuchos::RCP > solver); + virtual Teuchos::RCP > getSolver() const + { return solver_; } virtual void setObserver( Teuchos::RCP > obs = Teuchos::null); diff --git a/packages/tempus/src/Tempus_StepperIMEX_RK_decl.hpp b/packages/tempus/src/Tempus_StepperIMEX_RK_decl.hpp index d4e54769d449..1a74ace838b0 100644 --- a/packages/tempus/src/Tempus_StepperIMEX_RK_decl.hpp +++ b/packages/tempus/src/Tempus_StepperIMEX_RK_decl.hpp @@ -283,6 +283,8 @@ class StepperIMEX_RK : virtual public Tempus::StepperImplicit /// Set solver. virtual void setSolver( Teuchos::RCP > solver); + virtual Teuchos::RCP > getSolver() const + { return solver_; } virtual void setObserver( Teuchos::RCP > obs = Teuchos::null); diff --git a/packages/tempus/src/Tempus_StepperLeapfrog_decl.hpp b/packages/tempus/src/Tempus_StepperLeapfrog_decl.hpp index 30e0de508273..19137180a130 100644 --- a/packages/tempus/src/Tempus_StepperLeapfrog_decl.hpp +++ b/packages/tempus/src/Tempus_StepperLeapfrog_decl.hpp @@ -91,6 +91,8 @@ class StepperLeapfrog : virtual public Tempus::Stepper Teuchos::RCP solverPL=Teuchos::null); virtual void setSolver( Teuchos::RCP > solver); + virtual Teuchos::RCP > getSolver() const + { return Teuchos::null; } virtual void setObserver( Teuchos::RCP > obs = Teuchos::null); diff --git a/packages/tempus/src/Tempus_StepperNewmarkExplicitAForm_decl.hpp b/packages/tempus/src/Tempus_StepperNewmarkExplicitAForm_decl.hpp index 991fdaf2b7be..856726c2c428 100644 --- a/packages/tempus/src/Tempus_StepperNewmarkExplicitAForm_decl.hpp +++ b/packages/tempus/src/Tempus_StepperNewmarkExplicitAForm_decl.hpp @@ -43,6 +43,8 @@ class StepperNewmarkExplicitAForm : virtual public Tempus::Stepper Teuchos::RCP solverPL=Teuchos::null); virtual void setSolver( Teuchos::RCP > solver); + virtual Teuchos::RCP > getSolver() const + { return Teuchos::null; } /// Initialize during construction and after changing input parameters. virtual void initialize(){} diff --git a/packages/tempus/src/Tempus_StepperNewmarkImplicitAForm_decl.hpp b/packages/tempus/src/Tempus_StepperNewmarkImplicitAForm_decl.hpp index c03f54205f86..958c9a125e29 100644 --- a/packages/tempus/src/Tempus_StepperNewmarkImplicitAForm_decl.hpp +++ b/packages/tempus/src/Tempus_StepperNewmarkImplicitAForm_decl.hpp @@ -63,6 +63,8 @@ class StepperNewmarkImplicitAForm Teuchos::RCP solverPL=Teuchos::null); virtual void setSolver( Teuchos::RCP > solver); + virtual Teuchos::RCP > getSolver() const + { return solver_; } /// Initialize during construction and after changing input parameters. diff --git a/packages/tempus/src/Tempus_StepperNewmarkImplicitDForm_decl.hpp b/packages/tempus/src/Tempus_StepperNewmarkImplicitDForm_decl.hpp index 67e1f06ee9bd..f9bf94f3085b 100644 --- a/packages/tempus/src/Tempus_StepperNewmarkImplicitDForm_decl.hpp +++ b/packages/tempus/src/Tempus_StepperNewmarkImplicitDForm_decl.hpp @@ -62,6 +62,8 @@ class StepperNewmarkImplicitDForm : virtual public Tempus::StepperImplicit solverPL = Teuchos::null); virtual void setSolver(Teuchos::RCP> solver); + virtual Teuchos::RCP > getSolver() const + { return solver_; } /// Initialize during construction and after changing input parameters. virtual void diff --git a/packages/tempus/src/Tempus_StepperStaggeredForwardSensitivity.cpp b/packages/tempus/src/Tempus_StepperStaggeredForwardSensitivity.cpp new file mode 100644 index 000000000000..a10b4cc9ee83 --- /dev/null +++ b/packages/tempus/src/Tempus_StepperStaggeredForwardSensitivity.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_StepperStaggeredForwardSensitivity_decl.hpp" +#include "Tempus_StepperStaggeredForwardSensitivity_impl.hpp" + +namespace Tempus { + TEMPUS_INSTANTIATE_TEMPLATE_CLASS(StepperStaggeredForwardSensitivity) +} + +#endif diff --git a/packages/tempus/src/Tempus_StepperStaggeredForwardSensitivity_decl.hpp b/packages/tempus/src/Tempus_StepperStaggeredForwardSensitivity_decl.hpp new file mode 100644 index 000000000000..e2c4fd207580 --- /dev/null +++ b/packages/tempus/src/Tempus_StepperStaggeredForwardSensitivity_decl.hpp @@ -0,0 +1,143 @@ +// @HEADER +// **************************************************************************** +// Tempus: Copyright (2017) Sandia Corporation +// +// Distributed under BSD 3-clause license (See accompanying file Copyright.txt) +// **************************************************************************** +// @HEADER + +#ifndef Tempus_StepperStaggeredForwardSensitivity_decl_hpp +#define Tempus_StepperStaggeredForwardSensitivity_decl_hpp + +#include "Tempus_Stepper.hpp" +#include "Tempus_StaggeredForwardSensitivityModelEvaluator.hpp" + +namespace Tempus { + +/** \brief A stepper implementing staggered forward sensitivity analysis. + */ +/** + * It constructs two internal steppers, one for the state equations as usual + * and one for the sensitivity equations using + * Tempus::StaggeredForwardSensitivityModelEvaluator. It's implementation + * of takeStep() first takes a step using the state stepper, updates the + * sensitivity model evaluator with the compute state solution and time + * derivatives, and then takes a step using the sensitivity stepper. It + * optionally can reuse the state solver for the sensitivity equations as well. + */ +template +class StepperStaggeredForwardSensitivity : + virtual public Tempus::Stepper +{ +public: + + /// Constructor + /*! + * The first parameter list argument supplies supplies regular stepper + * options, while the second provides sensitivity specific options: + *
    + *
  • "Reuse State Linear Solver" (default: false) Whether to reuse the + * model's W matrix, solver, and preconditioner when solving the + * sensitivity equations. If they can be reused, substantial savings + * in compute time are possible. + *
  • "Force W Update" (default: false) When reusing the solver as above + * whether to force recomputation of W. This can be necessary when + * the solver overwrites it during the solve-phase (e.g., by a + * factorization). + *
  • "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 in the Tempus::CombinedForwardSensitivityModelEvaluator + * documentation. + *
  • "Sensitivity Parameter Index" (default: 0) Model evaluator + * parameter index for which sensitivities will be computed. + *
  • "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. + *
  • "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. + *
  • "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). + *
+ */ + StepperStaggeredForwardSensitivity( + const Teuchos::RCP >& appModel, + const Teuchos::RCP& pList = Teuchos::null, + const Teuchos::RCP& sens_pList = Teuchos::null); + + /// \name Basic stepper methods + //@{ + virtual void setModel( + const Teuchos::RCP >& appModel); + virtual void setNonConstModel( + const Teuchos::RCP >& appModel); + virtual Teuchos::RCP > getModel(); + + virtual void setSolver(std::string solverName); + virtual void setSolver( + Teuchos::RCP solverPL=Teuchos::null); + virtual void setSolver( + Teuchos::RCP > solver); + virtual Teuchos::RCP > getSolver() const + { return stateStepper_->getSolver(); } + + /// Initialize during construction and after changing input parameters. + virtual void initialize(); + + /// Take the specified timestep, dt, and return true if successful. + virtual void takeStep( + const Teuchos::RCP >& solutionHistory); + + /// Get a default (initial) StepperState + virtual Teuchos::RCP > + getDefaultStepperState(); + virtual Scalar getOrder() const {return stateStepper_->getOrder();} + virtual Scalar getOrderMin() const {return stateStepper_->getOrderMin();} + virtual Scalar getOrderMax() const {return stateStepper_->getOrderMax();} + //@} + + /// \name ParameterList methods + //@{ + void setParameterList(const Teuchos::RCP & pl); + Teuchos::RCP getNonconstParameterList(); + Teuchos::RCP unsetParameterList(); + Teuchos::RCP getValidParameters() const; + Teuchos::RCP getDefaultParameters() const; + //@} + + /// \name Overridden from Teuchos::Describable + //@{ + virtual std::string description() const; + virtual void describe(Teuchos::FancyOStream & out, + const Teuchos::EVerbosityLevel verbLevel) const; + //@} + + Teuchos::RCP > get_x_space() const; + +private: + + /// Default Constructor -- not allowed + StepperStaggeredForwardSensitivity(); + + void setParams(const Teuchos::RCP & pl, + const Teuchos::RCP & spl); + +private: + + Teuchos::RCP stepperPL_; + Teuchos::RCP sensPL_; + Teuchos::RCP > stateStepper_; + Teuchos::RCP > sensitivityStepper_; + Teuchos::RCP > fsa_model_; + Teuchos::RCP > stateSolutionHistory_; + Teuchos::RCP > sensSolutionHistory_; + bool reuse_solver_; + bool force_W_update_; +}; + +} // namespace Tempus + +#endif // Tempus_StepperStaggeredForwardSensitivity_decl_hpp diff --git a/packages/tempus/src/Tempus_StepperStaggeredForwardSensitivity_impl.hpp b/packages/tempus/src/Tempus_StepperStaggeredForwardSensitivity_impl.hpp new file mode 100644 index 000000000000..2ea156d16a09 --- /dev/null +++ b/packages/tempus/src/Tempus_StepperStaggeredForwardSensitivity_impl.hpp @@ -0,0 +1,425 @@ +// @HEADER +// **************************************************************************** +// Tempus: Copyright (2017) Sandia Corporation +// +// Distributed under BSD 3-clause license (See accompanying file Copyright.txt) +// **************************************************************************** +// @HEADER + +#ifndef Tempus_StepperStaggeredForwardSensitivity_impl_hpp +#define Tempus_StepperStaggeredForwardSensitivity_impl_hpp + +#include "Tempus_config.hpp" +#include "Tempus_StepperFactory.hpp" +#include "Tempus_StaggeredForwardSensitivityModelEvaluator.hpp" +#include "Teuchos_VerboseObjectParameterListHelpers.hpp" + +#include "Thyra_DefaultMultiVectorProductVectorSpace.hpp" +#include "Thyra_DefaultMultiVectorProductVector.hpp" +#include "NOX_Thyra.H" + +namespace Tempus { + +// Forward Declaration for recursive includes (this Stepper <--> StepperFactory) +template class StepperFactory; + +// StepperStaggeredForwardSensitivity definitions: +template +StepperStaggeredForwardSensitivity:: +StepperStaggeredForwardSensitivity( + const Teuchos::RCP >& appModel, + const Teuchos::RCP& pList, + const Teuchos::RCP& sens_pList) +{ + using Teuchos::RCP; + using Teuchos::ParameterList; + + // Set all the input parameters and call initialize + this->setParams(pList, sens_pList); + this->setModel(appModel); + this->initialize(); +} + + +template +void StepperStaggeredForwardSensitivity:: +setModel( + const Teuchos::RCP >& appModel) +{ + using Teuchos::RCP; + using Teuchos::rcp; + using Teuchos::ParameterList; + + // Create forward sensitivity model evaluator wrapper + Teuchos::RCP spl = Teuchos::parameterList(); + *spl = *sensPL_; + spl->remove("Reuse State Linear Solver"); + spl->remove("Force W Update"); + fsa_model_ = + rcp(new StaggeredForwardSensitivityModelEvaluator( + appModel, spl)); + + // Create state and sensitivity steppers + RCP > sf =Teuchos::rcp(new StepperFactory()); + if (stateStepper_ == Teuchos::null) + stateStepper_ = sf->createStepper(appModel, stepperPL_); + else + stateStepper_->setModel(appModel); + if (sensitivityStepper_ == Teuchos::null) + sensitivityStepper_ = sf->createStepper(fsa_model_, stepperPL_); + else + sensitivityStepper_->setModel(fsa_model_); +} + + +template +void StepperStaggeredForwardSensitivity:: +setNonConstModel( + const Teuchos::RCP >& appModel) +{ + this->setModel(appModel); +} + + +template +void StepperStaggeredForwardSensitivity:: +setSolver(std::string solverName) +{ + stateStepper_->setSolver(solverName); + sensitivityStepper_->setSolver(solverName); +} + +template +Teuchos::RCP > +StepperStaggeredForwardSensitivity:: +getModel() +{ + return stateStepper_->getModel(); +} + + +template +void StepperStaggeredForwardSensitivity:: +setSolver( + Teuchos::RCP solverPL) +{ + stateStepper_->setSolver(solverPL); + sensitivityStepper_->setSolver(solverPL); +} + + +template +void StepperStaggeredForwardSensitivity:: +setSolver( + Teuchos::RCP > solver) +{ + stateStepper_->setSolver(solver); + sensitivityStepper_->setSolver(solver); +} + + +template +void StepperStaggeredForwardSensitivity:: +initialize() +{ + this->setSolver(); +} + + +template +void StepperStaggeredForwardSensitivity:: +takeStep( + const Teuchos::RCP >& solutionHistory) +{ + using Teuchos::RCP; + using Teuchos::rcp; + using Teuchos::rcp_dynamic_cast; + using Thyra::VectorBase; + using Thyra::MultiVectorBase; + using Thyra::assign; + using Thyra::createMember; + using Thyra::multiVectorProductVector; + using Thyra::multiVectorProductVectorSpace; + typedef Thyra::DefaultMultiVectorProductVector DMVPV; + typedef Thyra::DefaultMultiVectorProductVectorSpace DMVPVS; + + // Initialize state, sensitivity solution histories if necessary. + // We only need to split the solution history into state and sensitivity + // components for the first step, otherwise the state and sensitivity + // histories are updated from the previous step. + if (stateSolutionHistory_ == Teuchos::null) { + RCP shPL = + solutionHistory->getNonconstParameterList(); + + // Get product X, XDot, XDotDot + RCP > state = solutionHistory->getCurrentState(); + RCP X, XDot, XDotDot; + X = rcp_dynamic_cast(state->getX()); + XDot = rcp_dynamic_cast(state->getXDot()); + if (state->getXDotDot() != Teuchos::null) + XDotDot = rcp_dynamic_cast(state->getXDotDot()); + + // Pull out state components (has to be non-const because of SolutionState + // constructor) + RCP > x, xdot, xdotdot; + x = X->getNonconstMultiVector()->col(0); + xdot = XDot->getNonconstMultiVector()->col(0); + if (XDotDot != Teuchos::null) + xdotdot = XDotDot->getNonconstMultiVector()->col(0); + + // Create state solution history + RCP > state_state = + rcp(new SolutionState(state->getMetaData()->clone(), + x, xdot, xdotdot, + state->getStepperState()->clone())); + stateSolutionHistory_ = rcp(new SolutionHistory(shPL)); + stateSolutionHistory_->addState(state_state); + + const int num_param = X->getMultiVector()->domain()->dim()-1; + const Teuchos::Range1D rng(1,num_param); + + // Pull out sensitivity components + RCP > dxdp, dxdotdp, dxdotdotdp; + dxdp = X->getNonconstMultiVector()->subView(rng); + dxdotdp = XDot->getNonconstMultiVector()->subView(rng); + if (XDotDot != Teuchos::null) + dxdotdotdp = XDotDot->getNonconstMultiVector()->subView(rng); + + // Create sensitivity product vectors + RCP dxdp_vec, dxdotdp_vec, dxdotdotdp_vec; + RCP prod_space = + multiVectorProductVectorSpace(X->getMultiVector()->range(), num_param); + dxdp_vec = multiVectorProductVector(prod_space, dxdp); + dxdotdp_vec = multiVectorProductVector(prod_space, dxdotdp); + if (XDotDot != Teuchos::null) + dxdotdotdp_vec = multiVectorProductVector(prod_space, dxdotdotdp); + + // Create sensitivity solution history + RCP > sens_state = + rcp(new SolutionState(state->getMetaData()->clone(), + dxdp_vec, dxdotdp_vec, dxdotdotdp_vec, + state->getStepperState()->clone())); + sensSolutionHistory_ = rcp(new SolutionHistory(shPL)); + sensSolutionHistory_->addState(sens_state); + } + + // Take step for state equations + stateSolutionHistory_->initWorkingState(); + stateSolutionHistory_->getWorkingState()->getMetaData()->copy( + solutionHistory->getWorkingState()->getMetaData()); + stateStepper_->takeStep(stateSolutionHistory_); + + // Get current state and set in sensitivity model evaluator + RCP > state = + stateSolutionHistory_->getWorkingState(); + RCP > x = state->getX(); + RCP > xdot = state->getXDot(); + RCP > xdotdot = state->getXDotDot(); + fsa_model_->setModelX(x); + fsa_model_->setModelXDot(xdot); + fsa_model_->setModelXDotDot(xdotdot); + + // Reuse state solver if requested + if (reuse_solver_ && stateStepper_->getSolver() != Teuchos::null) { + RCP nox_solver = + rcp_dynamic_cast(stateStepper_->getSolver()); + fsa_model_->setLO(nox_solver->get_nonconst_W_op(force_W_update_)); + fsa_model_->setPO(nox_solver->get_nonconst_prec_op()); + } + + // Take step in senstivity equations + sensSolutionHistory_->initWorkingState(); + sensSolutionHistory_->getWorkingState()->getMetaData()->copy( + solutionHistory->getWorkingState()->getMetaData()); + sensitivityStepper_->takeStep(sensSolutionHistory_); + + // Get current sensitivities + RCP > dxdp, dxdotdp, dxdotdotdp; + RCP > sens_state = + sensSolutionHistory_->getWorkingState(); + dxdp = rcp_dynamic_cast( + sens_state->getX())->getMultiVector(); + dxdotdp = rcp_dynamic_cast( + sens_state->getXDot())->getMultiVector(); + if (sens_state->getXDotDot() != Teuchos::null) + dxdotdotdp = rcp_dynamic_cast( + sens_state->getXDotDot())->getMultiVector(); + + const int num_param = dxdp->domain()->dim(); + const Teuchos::Range1D rng(1,num_param); + + // Form product solutions + RCP X, XDot, XDotDot; + RCP > prod_space = + multiVectorProductVectorSpace(dxdp->range(), num_param+1); + X = rcp_dynamic_cast(createMember(prod_space)); + assign(X->getNonconstMultiVector()->col(0).ptr(), *x); + assign(X->getNonconstMultiVector()->subView(rng).ptr(), *dxdp); + XDot = rcp_dynamic_cast(createMember(prod_space)); + assign(XDot->getNonconstMultiVector()->col(0).ptr(), *xdot); + assign(XDot->getNonconstMultiVector()->subView(rng).ptr(), *dxdotdp); + if (xdotdot != Teuchos::null) { + XDotDot = rcp_dynamic_cast(createMember(prod_space)); + assign(XDotDot->getNonconstMultiVector()->col(0).ptr(), *xdot); + assign(XDotDot->getNonconstMultiVector()->subView(rng).ptr(), *dxdotdotdp); + } + + // Add step to solution history + RCP > prod_state = solutionHistory->getWorkingState(); + assign(prod_state->getX().ptr(), *X); + assign(prod_state->getXDot().ptr(), *XDot); + if (XDotDot != Teuchos::null) + assign(prod_state->getXDotDot().ptr(), *XDotDot); + prod_state->setOrder(std::min(state->getOrder(), sens_state->getOrder())); + + // Determine whether step passed or failed + bool passed = true; + if (state->getStepperStatus() == Status::FAILED || + state->getSolutionStatus() == Status::FAILED || + sens_state->getStepperStatus() == Status::FAILED || + sens_state->getSolutionStatus() == Status::FAILED) + passed = false; + + if (passed) { + prod_state->getStepperState()->stepperStatus_ = Status::PASSED; + stateSolutionHistory_->promoteWorkingState(); + sensSolutionHistory_->promoteWorkingState(); + } + else + prod_state->getStepperState()->stepperStatus_ = Status::FAILED; +} + + +template +Teuchos::RCP > +StepperStaggeredForwardSensitivity:: +getDefaultStepperState() +{ + // ETP: Note, maybe this should be a special stepper state that combines + // states from both state and sensitivity steppers? + Teuchos::RCP > stepperState = + rcp(new StepperState(description())); + return stepperState; +} + + +template +std::string StepperStaggeredForwardSensitivity:: +description() const +{ + std::string name = "StepperStaggeredForwardSensitivity"; + return(name); +} + + +template +void StepperStaggeredForwardSensitivity:: +describe( + Teuchos::FancyOStream &out, + const Teuchos::EVerbosityLevel verbLevel) const +{ + out << description() << "::describe:" << std::endl; + stateStepper_->describe(out, verbLevel); + out << std::endl; + sensitivityStepper_->describe(out, verbLevel); +} + + +template +void StepperStaggeredForwardSensitivity:: +setParameterList( + Teuchos::RCP const& pList) +{ + if (pList == Teuchos::null) + stepperPL_ = this->getDefaultParameters(); + else + stepperPL_ = pList; + // Can not validate because of optional Parameters (e.g., Solver Name). + //stepperPL_->validateParametersAndSetDefaults(*this->getValidParameters()); +} + + +template +Teuchos::RCP +StepperStaggeredForwardSensitivity:: +getValidParameters() const +{ + return stateStepper_->getValidParameters(); +} + + +template +Teuchos::RCP +StepperStaggeredForwardSensitivity:: +getDefaultParameters() const +{ + return stateStepper_->getDefaultParameters(); +} + + +template +Teuchos::RCP +StepperStaggeredForwardSensitivity:: +getNonconstParameterList() +{ + return stepperPL_; +} + + +template +Teuchos::RCP +StepperStaggeredForwardSensitivity:: +unsetParameterList() +{ + Teuchos::RCP temp_plist = stepperPL_; + stepperPL_ = Teuchos::null; + return temp_plist; +} + + +template +void StepperStaggeredForwardSensitivity:: +setParams( + Teuchos::RCP const& pList, + Teuchos::RCP const& spList) +{ + if (pList == Teuchos::null) + stepperPL_ = this->getDefaultParameters(); + else + stepperPL_ = pList; + + if (spList == Teuchos::null) + sensPL_ = Teuchos::parameterList(); + else + sensPL_ = spList; + + reuse_solver_ = sensPL_->get("Reuse State Linear Solver", false); + force_W_update_ = sensPL_->get("Force W Update", true); + + // Can not validate because of optional Parameters (e.g., Solver Name). + //stepperPL_->validateParametersAndSetDefaults(*this->getValidParameters()); +} + + +template +Teuchos::RCP > +StepperStaggeredForwardSensitivity:: +get_x_space() const +{ + using Teuchos::RCP; + using Teuchos::rcp_dynamic_cast; + typedef Thyra::DefaultMultiVectorProductVectorSpace DMVPVS; + + RCP > x_space = + stateStepper_->getModel()->get_x_space(); + RCP dxdp_space = + rcp_dynamic_cast(sensitivityStepper_->getModel()->get_x_space()); + const int num_param = dxdp_space->numBlocks(); + RCP > prod_space = + multiVectorProductVectorSpace(x_space, num_param+1); + return prod_space; +} + +} // namespace Tempus + +#endif // Tempus_StepperStaggeredForwardSensitivity_impl_hpp diff --git a/packages/tempus/src/Thyra_MultiVectorLinearOp.hpp b/packages/tempus/src/Thyra_MultiVectorLinearOp.hpp new file mode 100644 index 000000000000..aeee19066819 --- /dev/null +++ b/packages/tempus/src/Thyra_MultiVectorLinearOp.hpp @@ -0,0 +1,359 @@ +// @HEADER +// **************************************************************************** +// Tempus: Copyright (2017) Sandia Corporation +// +// Distributed under BSD 3-clause license (See accompanying file Copyright.txt) +// **************************************************************************** +// @HEADER + +#ifndef Thyra_MultiVectorLinearOp_hpp +#define Thyra_MultiVectorLinearOp_hpp + +#include "Thyra_RowStatLinearOpBase.hpp" +#include "Thyra_ScaledLinearOpBase.hpp" +#include "Thyra_DefaultMultiVectorProductVectorSpace.hpp" +#include "Teuchos_ConstNonconstObjectContainer.hpp" + +#include "Thyra_DefaultMultiVectorProductVectorSpace.hpp" +#include "Thyra_AssertOp.hpp" +#include "Teuchos_dyn_cast.hpp" + +namespace Thyra { + +/** + * \brief Implicit concrete LinearOpBase subclass that + * takes a flattended out multi-vector and performs a multi-RHS apply with it. + */ +template +class MultiVectorLinearOp : + virtual public RowStatLinearOpBase, + virtual public ScaledLinearOpBase +{ +public: + + /** @name Constructors/initializers/accessors */ + //@{ + + /** \brief Construct to uninitialized. */ + MultiVectorLinearOp() {} + + /** \brief . */ + void nonconstInitialize( + const RCP > &op, + const RCP > &multiVecRange, + const RCP > &multiVecDomain + ) { + validateInitialize(op,multiVecRange,multiVecDomain); + op_ = op; + multiVecRange_ = multiVecRange; + multiVecDomain_ = multiVecDomain; + } + + /** \brief . */ + void initialize( + const RCP > &op, + const RCP > &multiVecRange, + const RCP > &multiVecDomain + ) { + validateInitialize(op,multiVecRange,multiVecDomain); + op_ = op; + multiVecRange_ = multiVecRange; + multiVecDomain_ = multiVecDomain; + } + + /** \brief . */ + RCP > + getNonconstLinearOp() { return op_.getNonconstObj(); } + + /** \brief . */ + RCP > + getLinearOp() const { return op_.getConstObj(); } + + /** \brief . */ + void uninitialize() { + op_.uninitialize(); + multiVecRange_ = Teuchos::null; + multiVecDomain_ = Teuchos::null; + } + + //@} + + /** @name Overridden from LinearOpBase */ + //@{ + + /** \brief . */ + RCP< const VectorSpaceBase > range() const { return multiVecRange_; } + + /** \brief . */ + RCP< const VectorSpaceBase > domain() const { return multiVecDomain_; } + + /** \brief . */ + RCP > clone() const { return Teuchos::null; // ToDo: Implement if needed ??? + } + //@} + +protected: + + /** @name Overridden from LinearOpBase */ + //@{ + /** \brief . */ + bool opSupportedImpl(EOpTransp M_trans) const { + return Thyra::opSupported(*op_.getConstObj(),M_trans); + } + + /** \brief . */ + void applyImpl( + const EOpTransp M_trans, + const MultiVectorBase &XX, + const Ptr > &YY, + const Scalar alpha, + const Scalar beta + ) const { + + using Teuchos::dyn_cast; + typedef DefaultMultiVectorProductVector MVPV; + + const Ordinal numCols = XX.domain()->dim(); + + for (Ordinal col_j = 0; col_j < numCols; ++col_j) { + + const RCP > x = XX.col(col_j); + const RCP > y = YY->col(col_j); + + RCP > + X = dyn_cast(*x).getMultiVector().assert_not_null(); + RCP > + Y = dyn_cast(*y).getNonconstMultiVector().assert_not_null(); + + Thyra::apply( *op_.getConstObj(), M_trans, *X, Y.ptr(), alpha, beta ); + } + + } + //@} + + /** @name Overridden from RowStatLinearOpBase */ + //@{ + + /** \brief Determine if a given row stat is supported. */ + bool rowStatIsSupportedImpl( + const RowStatLinearOpBaseUtils::ERowStat rowStat + ) const + { using Teuchos::rcp_dynamic_cast; + return rcp_dynamic_cast >(op_.getConstObj())->rowStatIsSupported(rowStat); } + + /** \brief Get some statistics about a supported row. + * + * \precondition this->rowStatIsSupported(rowStat)==true + */ + void getRowStatImpl( + const RowStatLinearOpBaseUtils::ERowStat rowStat, + const Ptr > &rowStatVec + ) const + { + TEUCHOS_ASSERT(this->rowStatIsSupported(rowStat)); + + // Compute the scaling vector from the underlying operator and assign + // it to each column of the scaling multivector. We only use the first + // column in the scaling below, but this makes the scaling vector + // consistent with a Kronecker-product operator + typedef DefaultMultiVectorProductVector MVPV; + using Teuchos::dyn_cast; + using Teuchos::rcp_dynamic_cast; + RCP > rowStatMultiVec = + dyn_cast(*rowStatVec).getNonconstMultiVector().assert_not_null(); + const Ordinal numCols = rowStatMultiVec->domain()->dim(); + if (numCols > 0) { + rcp_dynamic_cast >(op_.getConstObj())->getRowStat(rowStat, rowStatMultiVec->col(0).ptr()); + for (Ordinal col=1; colcol(0)), + rowStatMultiVec->col(col).ptr()); + } + } + } + + //@} + + /** @name Overridden from ScaledLinearOpBase */ + //@{ + + /** \brief . */ + virtual bool supportsScaleLeftImpl() const { + using Teuchos::rcp_dynamic_cast; + return rcp_dynamic_cast >(op_.getConstObj())->supportsScaleLeft(); + } + + /** \brief . */ + virtual bool supportsScaleRightImpl() const { + using Teuchos::rcp_dynamic_cast; + return rcp_dynamic_cast >(op_.getConstObj())->supportsScaleRight(); + } + + /** \brief . */ + virtual void scaleLeftImpl(const VectorBase &row_scaling) { + TEUCHOS_ASSERT(this->supportsScaleLeft()); + + // Use the first column in the scaling vector to scale the operator + typedef DefaultMultiVectorProductVector MVPV; + using Teuchos::dyn_cast; + using Teuchos::rcp_dynamic_cast; + RCP > row_scaling_mv = + dyn_cast(row_scaling).getMultiVector().assert_not_null(); + const Ordinal numCols = row_scaling_mv->domain()->dim(); + if (numCols > 0) { + rcp_dynamic_cast >(op_.getNonconstObj())->scaleLeft(*(row_scaling_mv->col(0))); + } + + //row_scaling.describe(std::cout, Teuchos::VERB_EXTREME); + } + + /** \brief . */ + virtual void scaleRightImpl(const VectorBase &col_scaling) { + TEUCHOS_ASSERT(this->supportsScaleRight()); + + // Use the first column in the scaling vector to scale the operator + typedef DefaultMultiVectorProductVector MVPV; + using Teuchos::dyn_cast; + using Teuchos::rcp_dynamic_cast; + RCP > col_scaling_mv = + dyn_cast(col_scaling).getMultiVector().assert_not_null(); + const Ordinal numCols = col_scaling_mv->domain()->dim(); + if (numCols > 0) { + rcp_dynamic_cast >(op_.getNonconstObj())->scaleRight(*(col_scaling_mv->col(0))); + } + } + + //@} + +private: + + // ////////////////////////////// + // Private types + + typedef Teuchos::ConstNonconstObjectContainer > CNOP; + + // ////////////////////////////// + // Private data members + + CNOP op_; + RCP > multiVecRange_; + RCP > multiVecDomain_; + + // ////////////////////////////// + // Private member functions + + static void validateInitialize( + const RCP > &op, + const RCP > &multiVecRange, + const RCP > &multiVecDomain + ) { +#ifdef TEUCHOS_DEBUG + TEUCHOS_TEST_FOR_EXCEPT(is_null(op)); + TEUCHOS_TEST_FOR_EXCEPT(is_null(multiVecRange)); + TEUCHOS_TEST_FOR_EXCEPT(is_null(multiVecDomain)); + TEUCHOS_TEST_FOR_EXCEPT( multiVecRange->numBlocks() != multiVecDomain->numBlocks() ); + if (op->range() != Teuchos::null) + THYRA_ASSERT_VEC_SPACES( + "MultiVectorLinearOp::initialize(op,multiVecRange,multiVecDomain)", + *op->range(), *multiVecRange->getBlock(0) ); + if (op->domain() != Teuchos::null) + THYRA_ASSERT_VEC_SPACES( + "MultiVectorLinearOp::initialize(op,multiVecRange,multiVecDomain)", + *op->domain(), *multiVecDomain->getBlock(0) ); +#endif +} + +}; + +/** \brief Nonmember constructor function. + * + * \relates MultiVectorLinearOp + */ +template +RCP > +multiVectorLinearOp() +{ + return Teuchos::rcp(new MultiVectorLinearOp()); +} + +/** \brief Nonmember constructor function. + * + * \relates MultiVectorLinearOp + */ +template +RCP > +nonconstMultiVectorLinearOp( + const RCP > &op, + const RCP > &multiVecRange, + const RCP > &multiVecDomain + ) +{ + RCP > + mvop = Teuchos::rcp(new MultiVectorLinearOp()); + mvop->nonconstInitialize(op,multiVecRange,multiVecDomain); + return mvop; +} + +/** \brief Nonmember constructor function. + * + * \relates MultiVectorLinearOp + */ +template +RCP > +nonconstMultiVectorLinearOp( + const RCP > &op, + const int num_blocks + ) +{ + RCP > + mvop = Teuchos::rcp(new MultiVectorLinearOp()); + RCP > mv_domain = + Thyra::multiVectorProductVectorSpace(op->domain(), num_blocks); + RCP > mv_range = + Thyra::multiVectorProductVectorSpace(op->range(), num_blocks); + mvop->nonconstInitialize(op,mv_range,mv_domain); + return mvop; +} + +/** \brief Nonmember constructor function. + * + * \relates MultiVectorLinearOp + */ +template +RCP > +multiVectorLinearOp( + const RCP > &op, + const RCP > &multiVecRange, + const RCP > &multiVecDomain + ) +{ + RCP > + mvop = Teuchos::rcp(new MultiVectorLinearOp()); + mvop->initialize(op,multiVecRange,multiVecDomain); + return mvop; +} + +/** \brief Nonmember constructor function. + * + * \relates MultiVectorLinearOp + */ +template +RCP > +multiVectorLinearOp( + const RCP > &op, + const int num_blocks + ) +{ + RCP > + mvop = Teuchos::rcp(new MultiVectorLinearOp()); + RCP > mv_domain = + Thyra::multiVectorProductVectorSpace(op->domain(), num_blocks); + RCP > mv_range = + Thyra::multiVectorProductVectorSpace(op->range(), num_blocks); + mvop->initialize(op,mv_range,mv_domain); + return mvop; +} + +} // end namespace Thyra + + +#endif diff --git a/packages/tempus/src/Thyra_MultiVectorLinearOpWithSolveFactory.hpp b/packages/tempus/src/Thyra_MultiVectorLinearOpWithSolveFactory.hpp new file mode 100644 index 000000000000..e989f95d1f3a --- /dev/null +++ b/packages/tempus/src/Thyra_MultiVectorLinearOpWithSolveFactory.hpp @@ -0,0 +1,668 @@ +// @HEADER +// **************************************************************************** +// Tempus: Copyright (2017) Sandia Corporation +// +// Distributed under BSD 3-clause license (See accompanying file Copyright.txt) +// **************************************************************************** +// @HEADER + +#ifndef Thyra_MultiVectorLinearOpWithSolveFactory_hpp +#define Thyra_MultiVectorLinearOpWithSolveFactory_hpp + +#include "Thyra_LinearOpWithSolveFactoryBase.hpp" +#include "Thyra_MultiVectorLinearOp.hpp" +#include "Thyra_MultiVectorPreconditioner.hpp" +#include "Thyra_MultiVectorPreconditionerFactory.hpp" +#include "Thyra_DefaultMultiVectorLinearOpWithSolve.hpp" +#include "Thyra_DefaultMultiVectorProductVectorSpace.hpp" +#include "Thyra_DefaultLinearOpSource.hpp" + +namespace Thyra { + +/** \brief Create a LinearOpWithSolveFactory for a flattened-out multi-vector. + */ +template +class MultiVectorLinearOpWithSolveFactory + : virtual public LinearOpWithSolveFactoryBase +{ +public: + + /** @name Overridden from Constructors/Initializers/Accessors */ + //@{ + + /** \brief Construct to uninitialized. */ + MultiVectorLinearOpWithSolveFactory() {} + + /** \brief Initialize given a single non-const LOWSFB object. + * + * \param lowsf [in,persisting] The LOWSFB object that will be used to + * create the LOWSB object for the diagonal blocks. + * + * Preconditions:
    + *
  • !is_null(lowsf) + *
+ * + */ + void nonconstInitialize( + const RCP > &lowsf, + const RCP > &multiVecRange, + const RCP > &multiVecDomain + ); + + + /** \brief Initialize given a single const LOWSFB object. + * + * \param lowsf [in,persisting] The LOWSFB object that will be used to + * create the LOWSB object for the diagonal blocks. + * + * Preconditions:
    + *
  • !is_null(lowsf) + *
+ * + */ + void initialize( + const RCP > &lowsf, + const RCP > &multiVecRange, + const RCP > &multiVecDomain + ); + + /** \brief . */ + RCP > getUnderlyingLOWSF(); + + /** \brief . */ + RCP > getUnderlyingLOWSF() const; + + //@} + + /** \name Overridden from Teuchos::Describable. */ + //@{ + + /** \brief . */ + std::string description() const; + + //@} + + /** @name Overridden from ParameterListAcceptor (simple forwarding functions) */ + //@{ + + /** \brief . */ + void setParameterList(RCP const& paramList); + /** \brief . */ + RCP getNonconstParameterList(); + /** \brief . */ + RCP unsetParameterList(); + /** \brief . */ + RCP getParameterList() const; + /** \brief . */ + RCP getValidParameters() const; + + //@} + + /** \name Overridden from LinearOpWithSolveFactoyBase */ + //@{ + + /** \brief returns false. */ + virtual bool acceptsPreconditionerFactory() const; + + /** \brief Throws exception. */ + virtual void setPreconditionerFactory( + const RCP > &precFactory, + const std::string &precFactoryName + ); + + /** \brief Returns null . */ + virtual RCP > + getPreconditionerFactory() const; + + /** \brief Throws exception. */ + virtual void unsetPreconditionerFactory( + RCP > *precFactory, + std::string *precFactoryName + ); + + /** \brief . */ + virtual bool isCompatible( + const LinearOpSourceBase &fwdOpSrc + ) const; + + /** \brief . */ + virtual RCP > createOp() const; + + /** \brief . */ + virtual void initializeOp( + const RCP > &fwdOpSrc, + LinearOpWithSolveBase *Op, + const ESupportSolveUse supportSolveUse + ) const; + + /** \brief . */ + virtual void initializeAndReuseOp( + const RCP > &fwdOpSrc, + LinearOpWithSolveBase *Op + ) const; + + /** \brief . */ + virtual void uninitializeOp( + LinearOpWithSolveBase *Op, + RCP > *fwdOpSrc, + RCP > *prec, + RCP > *approxFwdOpSrc, + ESupportSolveUse *supportSolveUse + ) const; + + /** \brief . */ + virtual bool supportsPreconditionerInputType( + const EPreconditionerInputType precOpType + ) const; + + /** \brief . */ + virtual void initializePreconditionedOp( + const RCP > &fwdOpSrc, + const RCP > &prec, + LinearOpWithSolveBase *Op, + const ESupportSolveUse supportSolveUse + ) const; + + /** \brief . */ + virtual void initializeApproxPreconditionedOp( + const RCP > &fwdOpSrc, + const RCP > &approxFwdOpSrc, + LinearOpWithSolveBase *Op, + const ESupportSolveUse supportSolveUse + ) const; + + //@} + +protected: + + /** \brief Overridden from Teuchos::VerboseObjectBase */ + //@{ + + /** \brief . */ + void informUpdatedVerbosityState() const; + + //@} + +private: + + typedef Teuchos::ConstNonconstObjectContainer > LOWSF_t; + + LOWSF_t lowsf_; + RCP > multiVecRange_; + RCP > multiVecDomain_; + +}; + +/** \brief Nonmember constructor. + * + * \releates MultiVectorLinearOpWithSolveFactory + */ +template +RCP > +nonconstMultiVectorLinearOpWithSolveFactory( + const RCP > &lowsf, + const RCP > &multiVecRange, + const RCP > &multiVecDomain + ) +{ + RCP > mvlowsf = + Teuchos::rcp(new MultiVectorLinearOpWithSolveFactory); + mvlowsf->nonconstInitialize(lowsf,multiVecRange,multiVecDomain); + return mvlowsf; +} + +/** \brief Nonmember constructor. + * + * \releates MultiVectorLinearOpWithSolveFactory + */ +template +RCP > +nonconstMultiVectorLinearOpWithSolveFactory( + const RCP > &lowsf, + const int num_blocks + ) +{ + RCP< LinearOpWithSolveBase > op = lowsf->createOp(); + RCP > mv_domain = + Thyra::multiVectorProductVectorSpace(op->domain(), num_blocks); + RCP > mv_range = + Thyra::multiVectorProductVectorSpace(op->range(), num_blocks); + return nonconstMultiVectorLinearOpWithSolveFactory(lowsf, mv_range, mv_domain); +} + +/** \brief Nonmember constructor. + * + * \releates MultiVectorLinearOpWithSolveFactory + */ +template +RCP > +multiVectorLinearOpWithSolveFactory( + const RCP > &lowsf, + const RCP > &multiVecRange, + const RCP > &multiVecDomain + ) +{ + RCP > mvlowsf = + Teuchos::rcp(new MultiVectorLinearOpWithSolveFactory); + mvlowsf->initialize(lowsf,multiVecRange,multiVecDomain); + return mvlowsf; +} + +/** \brief Nonmember constructor. + * + * \releates MultiVectorLinearOpWithSolveFactory + */ +template +RCP > +multiVectorLinearOpWithSolveFactory( + const RCP > &lowsf, + const int num_blocks + ) +{ + RCP< LinearOpWithSolveBase > op = lowsf->createOp(); + RCP > mv_domain = + Thyra::multiVectorProductVectorSpace(op->domain(), num_blocks); + RCP > mv_range = + Thyra::multiVectorProductVectorSpace(op->range(), num_blocks); + return multiVectorLinearOpWithSolveFactory(lowsf, mv_range, mv_domain); +} + +// Overridden from Constructors/Initializers/Accessors + +template +void +MultiVectorLinearOpWithSolveFactory:: +nonconstInitialize( + const RCP > &lowsf, + const RCP > &multiVecRange, + const RCP > &multiVecDomain + ) +{ +#ifdef TEUCHOS_DEBUG + TEUCHOS_TEST_FOR_EXCEPT(is_null(lowsf)); +#endif + lowsf_.initialize(lowsf); + multiVecRange_ = multiVecRange; + multiVecDomain_ = multiVecDomain; +} + +template +void +MultiVectorLinearOpWithSolveFactory:: +initialize( + const RCP > &lowsf, + const RCP > &multiVecRange, + const RCP > &multiVecDomain + ) +{ +#ifdef TEUCHOS_DEBUG + TEUCHOS_TEST_FOR_EXCEPT(is_null(lowsf)); +#endif + lowsf_.initialize(lowsf); + multiVecRange_ = multiVecRange; + multiVecDomain_ = multiVecDomain; +} + +template +RCP > +MultiVectorLinearOpWithSolveFactory:: +getUnderlyingLOWSF() +{ + return lowsf_.getNonconstObj(); +} + +template +RCP > +MultiVectorLinearOpWithSolveFactory:: +getUnderlyingLOWSF() const +{ + return lowsf_.getConstObj(); +} + +// Overridden from Teuchos::Describable + +template +std::string +MultiVectorLinearOpWithSolveFactory:: +description() const +{ + std::ostringstream oss; + oss << this->Teuchos::Describable::description() + << "{" + << "lowsf="; + if (!is_null(lowsf_.getConstObj())) + oss << lowsf_.getConstObj()->description(); + else + oss << "NULL"; + oss << "}"; + return oss.str(); +} + +// Overridden from ParameterListAcceptor + +template +void +MultiVectorLinearOpWithSolveFactory:: +setParameterList( + RCP const& paramList + ) +{ + lowsf_.getNonconstObj()->setParameterList(paramList); +} + +template +RCP +MultiVectorLinearOpWithSolveFactory:: +getNonconstParameterList() +{ + return lowsf_.getNonconstObj()->getNonconstParameterList(); +} + +template +RCP +MultiVectorLinearOpWithSolveFactory:: +unsetParameterList() +{ + return lowsf_.getNonconstObj()->unsetParameterList(); +} + +template +RCP +MultiVectorLinearOpWithSolveFactory:: +getParameterList() const +{ + return lowsf_.getConstObj()->getParameterList(); +} + +template +RCP +MultiVectorLinearOpWithSolveFactory:: +getValidParameters() const +{ + return lowsf_.getConstObj()->getValidParameters(); +} + +// Overridden from LinearOpWithSolveFactoyBase + +template +bool +MultiVectorLinearOpWithSolveFactory:: +acceptsPreconditionerFactory() const +{ + return lowsf_.getConstObj()->acceptsPreconditionerFactory(); +} + +template +void +MultiVectorLinearOpWithSolveFactory:: +setPreconditionerFactory( + const RCP > &precFactory, + const std::string &precFactoryName + ) +{ + typedef MultiVectorPreconditionerFactory MVPF; + RCP mvpf = Teuchos::rcp_dynamic_cast(precFactory); + lowsf_.getNonconstObj()->setPreconditionerFactory( + mvpf->getNonconstPreconditionerFactory(), + precFactoryName); +} + +template +RCP > +MultiVectorLinearOpWithSolveFactory:: +getPreconditionerFactory() const +{ + RCP > prec_fac = + lowsf_.getConstObj()->getPreconditionerFactory(); + if (prec_fac == Teuchos::null) + return Teuchos::null; + else + return nonconstMultiVectorPreconditionerFactory( + prec_fac, + multiVecRange_, multiVecDomain_); +} + +template +void MultiVectorLinearOpWithSolveFactory:: +unsetPreconditionerFactory( + RCP > *precFactory, + std::string *precFactoryName + ) +{ + RCP > inner_precFactory; + lowsf_.getNonconstObj()->unsetPreconditionerFactory( + precFactory ? &inner_precFactory : NULL, + precFactoryName); + if (precFactory) + *precFactory = nonconstMultiVectorPreconditionerFactory(inner_precFactory, + multiVecRange_, + multiVecDomain_); +} + +template +bool +MultiVectorLinearOpWithSolveFactory:: +isCompatible( + const LinearOpSourceBase &fwdOpSrc + ) const +{ + typedef MultiVectorLinearOp MVLO; + RCP mvlo = + Teuchos::rcp_dynamic_cast(fwdOpSrc.getOp().assert_not_null()); + RCP > inner_fwdOpSrc = + defaultLinearOpSource(mvlo->getLinearOp()); + return lowsf_.getConstObj()->isCompatible(*inner_fwdOpSrc); +} + +template +RCP > +MultiVectorLinearOpWithSolveFactory:: +createOp() const +{ + return nonconstMultiVectorLinearOpWithSolve( + lowsf_.getConstObj()->createOp(), + multiVecRange_, + multiVecDomain_); +} + +template +void +MultiVectorLinearOpWithSolveFactory:: +initializeOp( + const RCP > &fwdOpSrc, + LinearOpWithSolveBase *Op, + const ESupportSolveUse supportSolveUse + ) const +{ + using Teuchos::dyn_cast; + using Teuchos::rcp_dynamic_cast; + +#ifdef TEUCHOS_DEBUG + TEUCHOS_TEST_FOR_EXCEPT(0==Op); +#endif + + // Set the verbosity settings for the wrapped LOWSF object! + lowsf_.getConstObj()->setOStream(this->getOStream()); + lowsf_.getConstObj()->setVerbLevel(this->getVerbLevel()); + + typedef MultiVectorLinearOp MVLO; + typedef DefaultMultiVectorLinearOpWithSolve MVLOWS; + const RCP mvlo = + rcp_dynamic_cast(fwdOpSrc->getOp().assert_not_null()); + MVLOWS &mvlows = dyn_cast(*Op); + + lowsf_.getConstObj()->initializeOp( + defaultLinearOpSource(mvlo->getLinearOp()), + mvlows.getNonconstLinearOpWithSolve().get(), + supportSolveUse); +} + +template +void +MultiVectorLinearOpWithSolveFactory:: +initializeAndReuseOp( + const RCP > &fwdOpSrc, + LinearOpWithSolveBase *Op + ) const +{ + using Teuchos::dyn_cast; + using Teuchos::rcp_dynamic_cast; + +#ifdef TEUCHOS_DEBUG + TEUCHOS_TEST_FOR_EXCEPT(0==Op); +#endif + + // Set the verbosity settings for the wrapped LOWSF object! + lowsf_.getConstObj()->setOStream(this->getOStream()); + lowsf_.getConstObj()->setVerbLevel(this->getVerbLevel()); + + typedef MultiVectorLinearOp MVLO; + typedef DefaultMultiVectorLinearOpWithSolve MVLOWS; + const RCP mvlo = + rcp_dynamic_cast(fwdOpSrc->getOp().assert_not_null()); + MVLOWS &mvlows = dyn_cast(*Op); + + lowsf_.getConstObj()->initializeAndReuseOp( + defaultLinearOpSource(mvlo->getLinearOp()), + mvlows.getNonconstLinearOpWithSolve().get()); +} + +template +void +MultiVectorLinearOpWithSolveFactory:: +uninitializeOp( + LinearOpWithSolveBase *Op, + RCP > *fwdOpSrc, + RCP > *prec, + RCP > *approxFwdOpSrc, + ESupportSolveUse *supportSolveUse + ) const +{ + using Teuchos::dyn_cast; + +#ifdef TEUCHOS_DEBUG + TEUCHOS_TEST_FOR_EXCEPT(0==Op); +#endif + typedef DefaultMultiVectorLinearOpWithSolve MVLOWS; + MVLOWS &mvlowsOp = dyn_cast(*Op); + RCP > inner_fwdOpSrc; + RCP > inner_prec; + RCP > inner_approxFwdOpSrc; + lowsf_.getConstObj()->uninitializeOp( + mvlowsOp.getNonconstLinearOpWithSolve().get(), + fwdOpSrc ? &inner_fwdOpSrc : NULL, + prec ? &inner_prec : NULL, + approxFwdOpSrc ? &inner_approxFwdOpSrc : NULL, + supportSolveUse); + if (fwdOpSrc) + *fwdOpSrc = + defaultLinearOpSource(multiVectorLinearOp(inner_fwdOpSrc->getOp(), + multiVecRange_, + multiVecDomain_)); + if (prec) + *prec = multiVectorPreconditioner(inner_prec, + multiVecRange_, + multiVecDomain_); + if (fwdOpSrc) + *approxFwdOpSrc = + defaultLinearOpSource(multiVectorLinearOp(inner_approxFwdOpSrc->getOp(), + multiVecRange_, + multiVecDomain_)); +} + +template +bool +MultiVectorLinearOpWithSolveFactory:: +supportsPreconditionerInputType( + const EPreconditionerInputType precOpType + ) const +{ + return lowsf_.getConstObj()->supportsPreconditionerInputType(precOpType); +} + +template +void +MultiVectorLinearOpWithSolveFactory:: +initializePreconditionedOp( + const RCP > &fwdOpSrc, + const RCP > &prec, + LinearOpWithSolveBase *Op, + const ESupportSolveUse supportSolveUse + ) const +{ + using Teuchos::dyn_cast; + using Teuchos::rcp_dynamic_cast; + +#ifdef TEUCHOS_DEBUG + TEUCHOS_TEST_FOR_EXCEPT(0==Op); +#endif + + // Set the verbosity settings for the wrapped LOWSF object! + lowsf_.getConstObj()->setOStream(this->getOStream()); + lowsf_.getConstObj()->setVerbLevel(this->getVerbLevel()); + + typedef MultiVectorLinearOp MVLO; + typedef MultiVectorPreconditioner MVP; + typedef DefaultMultiVectorLinearOpWithSolve MVLOWS; + const RCP mvlo = + rcp_dynamic_cast(fwdOpSrc->getOp().assert_not_null()); + const RCP mvp = rcp_dynamic_cast(prec); + MVLOWS &mvlows = dyn_cast(*Op); + + lowsf_.getConstObj()->initializePreconditionedOp( + defaultLinearOpSource(mvlo->getLinearOp()), + mvp->getPreconditioner(), + mvlows.getNonconstLinearOpWithSolve().get(), + supportSolveUse); +} + +template +void +MultiVectorLinearOpWithSolveFactory:: +initializeApproxPreconditionedOp( + const RCP > &fwdOpSrc, + const RCP > &approxFwdOpSrc, + LinearOpWithSolveBase *Op, + const ESupportSolveUse supportSolveUse + ) const +{ + using Teuchos::dyn_cast; + using Teuchos::rcp_dynamic_cast; + +#ifdef TEUCHOS_DEBUG + TEUCHOS_TEST_FOR_EXCEPT(0==Op); +#endif + + // Set the verbosity settings for the wrapped LOWSF object! + lowsf_.getConstObj()->setOStream(this->getOStream()); + lowsf_.getConstObj()->setVerbLevel(this->getVerbLevel()); + + typedef MultiVectorLinearOp MVLO; + typedef DefaultMultiVectorLinearOpWithSolve MVLOWS; + const RCP mvlo = + rcp_dynamic_cast(fwdOpSrc->getOp().assert_not_null()); + const RCP amvlo = + rcp_dynamic_cast(approxFwdOpSrc->getOp().assert_not_null()); + MVLOWS &mvlows = dyn_cast(*Op); + + lowsf_.getConstObj()->initializeApproxPreconditionedOp( + defaultLinearOpSource(mvlo->getLinearOp()), + defaultLinearOpSource(amvlo->getLinearOp()), + mvlows.getNonconstLinearOpWithSolve().get(), + supportSolveUse); +} + +// protected + +template +void +MultiVectorLinearOpWithSolveFactory:: +informUpdatedVerbosityState() const +{ + lowsf_.getConstObj()->setVerbLevel(this->getVerbLevel()); + lowsf_.getConstObj()->setOStream(this->getOStream()); +} + +} // namespace Thyra + +#endif diff --git a/packages/tempus/src/Thyra_MultiVectorPreconditioner.hpp b/packages/tempus/src/Thyra_MultiVectorPreconditioner.hpp new file mode 100644 index 000000000000..82d01e24a6dd --- /dev/null +++ b/packages/tempus/src/Thyra_MultiVectorPreconditioner.hpp @@ -0,0 +1,213 @@ +// @HEADER +// **************************************************************************** +// Tempus: Copyright (2017) Sandia Corporation +// +// Distributed under BSD 3-clause license (See accompanying file Copyright.txt) +// **************************************************************************** +// @HEADER + +#ifndef Thyra_MultiVectorPreconditioner_hpp +#define Thyra_MultiVectorPreconditioner_hpp + +#include "Thyra_PreconditionerBase.hpp" +#include "Teuchos_ConstNonconstObjectContainer.hpp" +#include "Thyra_MultiVectorLinearOp.hpp" +#include "Thyra_DefaultMultiVectorProductVectorSpace.hpp" + +namespace Thyra { + +/** \brief Concrete PreconditionerBase subclass that + * wraps a preconditioner operator in MultiVectorLinearOp. + */ +template +class MultiVectorPreconditioner : virtual public PreconditionerBase +{ +public: + + /** @name Constructors/initializers/accessors */ + //@{ + + /** \brief Construct to uninitialized. */ + MultiVectorPreconditioner() {} + + /** \brief . */ + void nonconstInitialize( + const RCP > &prec, + const RCP > &multiVecRange, + const RCP > &multiVecDomain + ) { + validateInitialize(prec,multiVecRange,multiVecDomain); + prec_ = prec; + multiVecRange_ = multiVecRange; + multiVecDomain_ = multiVecDomain; + } + + /** \brief . */ + void initialize( + const RCP > &prec, + const RCP > &multiVecRange, + const RCP > &multiVecDomain) { + validateInitialize(prec,multiVecRange,multiVecDomain); + prec_ = prec; + multiVecRange_ = multiVecRange; + multiVecDomain_ = multiVecDomain; + } + + /** \brief . */ + RCP > + getNonconstPreconditioner() { return prec_.getNonconstObj(); } + + /** \brief . */ + RCP > + getPreconditioner() const { return prec_.getConstObj(); } + + /** \brief . */ + void uninitialize() { + prec_.uninitialize(); + multiVecRange_ = Teuchos::null; + multiVecDomain_ = Teuchos::null; + } + + //@} + + /** @name Overridden from PreconditionerBase */ + //@{ + + /** \brief . */ + bool isLeftPrecOpConst() const + { return prec_.getConstObj()->isLeftPrecOpConst(); } + + /** \brief . */ + Teuchos::RCP > getNonconstLeftPrecOp() + { return nonconstMultiVectorLinearOp( + prec_.getNonconstObj()->getNonconstLeftPrecOp(), + multiVecRange_, + multiVecDomain_); } + + /** \brief . */ + Teuchos::RCP > getLeftPrecOp() const + { return multiVectorLinearOp( + prec_.getConstObj()->getLeftPrecOp(), + multiVecRange_, + multiVecDomain_); } + + /** \brief . */ + bool isRightPrecOpConst() const + { return prec_.getConstObj()->isRightPrecOpConst(); } + + /** \brief . */ + Teuchos::RCP > getNonconstRightPrecOp() + { return nonconstMultiVectorLinearOp( + prec_.getNonconstObj()->getNonconstRightPrecOp(), + multiVecRange_, + multiVecDomain_); } + + /** \brief . */ + Teuchos::RCP > getRightPrecOp() const + { return multiVectorLinearOp( + prec_.getConstObj()->getRightPrecOp(), + multiVecRange_, + multiVecDomain_); } + + /** \brief . */ + bool isUnspecifiedPrecOpConst() const + { return prec_.getConstObj()->isUnspecifiedPrecOpConst(); } + + /** \brief . */ + Teuchos::RCP > getNonconstUnspecifiedPrecOp() + { return nonconstMultiVectorLinearOp( + prec_.getNonconstObj()->getNonconstUnspecifiedPrecOp(), + multiVecRange_, + multiVecDomain_); } + + /** \brief . */ + Teuchos::RCP > getUnspecifiedPrecOp() const + { return multiVectorLinearOp( + prec_.getNonconstObj()->getUnspecifiedPrecOp(), + multiVecRange_, + multiVecDomain_); } + + //@} + +private: + + // ////////////////////////////// + // Private types + + typedef Teuchos::ConstNonconstObjectContainer > CNPB; + + // ////////////////////////////// + // Private data members + + CNPB prec_; + RCP > multiVecRange_; + RCP > multiVecDomain_; + + // ////////////////////////////// + // Private member functions + + static void validateInitialize( + const RCP > &prec, + const RCP > &multiVecRange, + const RCP > &multiVecDomain + ) { +#ifdef TEUCHOS_DEBUG + TEUCHOS_TEST_FOR_EXCEPT(is_null(prec)); + TEUCHOS_TEST_FOR_EXCEPT(is_null(multiVecRange)); + TEUCHOS_TEST_FOR_EXCEPT(is_null(multiVecDomain)); + TEUCHOS_TEST_FOR_EXCEPT( multiVecRange->numBlocks() != multiVecDomain->numBlocks() ); +#endif + } + +}; + +/** \brief Nonmember constructor function. + * + * \relates MultiVectorPreconditioner + */ +template +RCP > +multiVectorPreconditioner() +{ + return Teuchos::rcp(new MultiVectorPreconditioner()); +} + +/** \brief Nonmember constructor function. + * + * \relates MultiVectorPreconditioner + */ +template +RCP > +nonconstMultiVectorPreconditioner( + const RCP > &prec, + const RCP > &multiVecRange, + const RCP > &multiVecDomain + ) +{ + RCP > + mvprec = Teuchos::rcp(new MultiVectorPreconditioner()); + mvprec->nonconstInitialize(prec,multiVecRange,multiVecDomain); + return mvprec; +} + +/** \brief Nonmember constructor function. + * + * \relates MultiVectorPreconditioner + */ +template +RCP > +multiVectorPreconditioner( + const RCP > &prec, + const RCP > &multiVecRange, + const RCP > &multiVecDomain + ) +{ + RCP > + mvprec = Teuchos::rcp(new MultiVectorPreconditioner()); + mvprec->initialize(prec,multiVecRange,multiVecDomain); + return mvprec; +} + +} // end namespace Thyra + +#endif diff --git a/packages/tempus/src/Thyra_MultiVectorPreconditionerFactory.hpp b/packages/tempus/src/Thyra_MultiVectorPreconditionerFactory.hpp new file mode 100644 index 000000000000..e1ffe7d91374 --- /dev/null +++ b/packages/tempus/src/Thyra_MultiVectorPreconditionerFactory.hpp @@ -0,0 +1,274 @@ +// @HEADER +// **************************************************************************** +// Tempus: Copyright (2017) Sandia Corporation +// +// Distributed under BSD 3-clause license (See accompanying file Copyright.txt) +// **************************************************************************** +// @HEADER + +#ifndef Thyra_MultiVectorPreconditionerFactory_hpp +#define Thyra_MultiVectorPreconditionerFactory_hpp + +#include "Thyra_PreconditionerFactoryBase.hpp" +#include "Teuchos_ConstNonconstObjectContainer.hpp" +#include "Thyra_MultiVectorLinearOp.hpp" +#include "Thyra_MultiVectorPreconditioner.hpp" +#include "Thyra_DefaultMultiVectorProductVectorSpace.hpp" + +namespace Thyra { + +/** \brief Concrete PreconditionerFactoryBase subclass that + * wraps a preconditioner in MultiVectorPreconditioner. + */ +template +class MultiVectorPreconditionerFactory + : virtual public PreconditionerFactoryBase +{ +public: + + /** @name Constructors/initializers/accessors */ + //@{ + + /** \brief Construct to uninitialized. */ + MultiVectorPreconditionerFactory() {} + + /** \brief . */ + void nonconstInitialize( + const RCP > &prec_fac, + const RCP > &multiVecRange, + const RCP > &multiVecDomain + ) { + validateInitialize(prec_fac,multiVecRange,multiVecDomain); + prec_fac_ = prec_fac; + multiVecRange_ = multiVecRange; + multiVecDomain_ = multiVecDomain; + } + + /** \brief . */ + void initialize( + const RCP > &prec_fac, + const RCP > &multiVecRange, + const RCP > &multiVecDomain) { + validateInitialize(prec_fac,multiVecRange,multiVecDomain); + prec_fac_ = prec_fac; + multiVecRange_ = multiVecRange; + multiVecDomain_ = multiVecDomain; + } + + /** \brief . */ + RCP > + getNonconstPreconditionerFactory() { return prec_fac_.getNonconstObj(); } + + /** \brief . */ + RCP > + getPreconditionerFactory() const { return prec_fac_.getConstObj(); } + + /** \brief . */ + void uninitialize() { + prec_fac_.uninitialize(); + multiVecRange_ = Teuchos::null; + multiVecDomain_ = Teuchos::null; + } + + /** \name Overridden from Teuchos::Describable. */ + //@{ + + /** \brief . */ + std::string description() const + { + std::ostringstream oss; + oss << this->Teuchos::Describable::description() + << "{" + << "prec_fac="; + if (!is_null(prec_fac_.getConstObj())) + oss << prec_fac_.getConstObj()->description(); + else + oss << "NULL"; + oss << "}"; + return oss.str(); + } + + //@} + + /** @name Overridden from ParameterListAcceptor (simple forwarding functions) */ + //@{ + + /** \brief . */ + void setParameterList(RCP const& paramList) + { + prec_fac_.getNonconstObj()->setParameterList(paramList); + } + + /** \brief . */ + RCP getNonconstParameterList() + { + return prec_fac_.getNonconstObj()->getNonconstParameterList(); + } + + /** \brief . */ + RCP unsetParameterList() + { + return prec_fac_.getNonconstObj()->unsetParameterList(); + } + + /** \brief . */ + RCP getParameterList() const + { + return prec_fac_.getConstObj()->getParameterList(); + } + + /** \brief . */ + RCP getValidParameters() const + { + return prec_fac_.getConstObj()->getValidParameters(); + } + + //@} + + //@} + + /** @name Overridden from PreconditionerFactoryBase */ + //@{ + + /** \brief . */ + bool isCompatible(const LinearOpSourceBase &fwdOpSrc) const + { return prec_fac_.getConstObj()->isCompatible(fwdOpSrc); } + + /** \brief . */ + RCP > createPrec() const + { return nonconstMultiVectorPreconditioner( + prec_fac_.getConstObj()->createPrec(), + multiVecRange_, + multiVecDomain_); } + + /** \brief . */ + void initializePrec( + const RCP > &fwdOpSrc, + PreconditionerBase *precOp, + const ESupportSolveUse supportSolveUse = SUPPORT_SOLVE_UNSPECIFIED + ) const + { + using Teuchos::dyn_cast; + using Teuchos::rcp_dynamic_cast; + + typedef MultiVectorLinearOp MVLO; + typedef MultiVectorPreconditioner MVP; + const RCP mvlo = + rcp_dynamic_cast(fwdOpSrc->getOp().assert_not_null()); + MVP &mvp = dyn_cast(*precOp); + prec_fac_.getConstObj()->initializePrec( + defaultLinearOpSource(mvlo->getLinearOp()), + mvp.getNonconstPreconditioner().get(), + supportSolveUse); + } + + /** \brief . */ + void uninitializePrec( + PreconditionerBase *precOp, + RCP > *fwdOpSrc = NULL, + ESupportSolveUse *supportSolveUse = NULL + ) const + { + using Teuchos::dyn_cast; + +#ifdef TEUCHOS_DEBUG + TEUCHOS_TEST_FOR_EXCEPT(0==precOp); +#endif + typedef MultiVectorPreconditioner MVP; + MVP &mvp = dyn_cast(*precOp); + RCP > inner_fwdOpSrc; + prec_fac_.getConstObj()->uninitializePrec( + mvp.getNonconstPreconditioner().get(), + fwdOpSrc ? &inner_fwdOpSrc : NULL, + supportSolveUse); + if (fwdOpSrc) + *fwdOpSrc = + defaultLinearOpSource(multiVectorLinearOp(inner_fwdOpSrc->getOp(), + multiVecRange_, + multiVecDomain_)); + } + + //@} + +private: + + // ////////////////////////////// + // Private types + + typedef Teuchos::ConstNonconstObjectContainer > CNPFB; + + // ////////////////////////////// + // Private data members + + CNPFB prec_fac_; + RCP > multiVecRange_; + RCP > multiVecDomain_; + + // ////////////////////////////// + // Private member functions + + static void validateInitialize( + const RCP > &prec_fac, + const RCP > &multiVecRange, + const RCP > &multiVecDomain + ) { +#ifdef TEUCHOS_DEBUG + TEUCHOS_TEST_FOR_EXCEPT(is_null(prec_fac)); + TEUCHOS_TEST_FOR_EXCEPT(is_null(multiVecRange)); + TEUCHOS_TEST_FOR_EXCEPT(is_null(multiVecDomain)); + TEUCHOS_TEST_FOR_EXCEPT( multiVecRange->numBlocks() != multiVecDomain->numBlocks() ); +#endif + } + +}; + +/** \brief Nonmember constructor function. + * + * \relates MultiVectorPreconditionerFactory + */ +template +RCP > +multiVectorPreconditionerFactory() +{ + return Teuchos::rcp(new MultiVectorPreconditionerFactory()); +} + +/** \brief Nonmember constructor function. + * + * \relates MultiVectorPreconditionerFactory + */ +template +RCP > +nonconstMultiVectorPreconditionerFactory( + const RCP > &prec_fac, + const RCP > &multiVecRange, + const RCP > &multiVecDomain + ) +{ + RCP > + mvfac = Teuchos::rcp(new MultiVectorPreconditionerFactory()); + mvfac->nonconstInitialize(prec_fac,multiVecRange,multiVecDomain); + return mvfac; +} + +/** \brief Nonmember constructor function. + * + * \relates MultiVectorPreconditionerFactory + */ +template +RCP > +multiVectorPreconditionerFactory( + const RCP > &prec_fac, + const RCP > &multiVecRange, + const RCP > &multiVecDomain + ) +{ + RCP > + mvfac = Teuchos::rcp(new MultiVectorPreconditionerFactory()); + mvfac->initialize(prec_fac,multiVecRange,multiVecDomain); + return mvfac; +} + +} // end namespace Thyra + +#endif diff --git a/packages/tempus/src/Thyra_ReuseLinearOpWithSolveFactory.hpp b/packages/tempus/src/Thyra_ReuseLinearOpWithSolveFactory.hpp new file mode 100644 index 000000000000..997e323956dd --- /dev/null +++ b/packages/tempus/src/Thyra_ReuseLinearOpWithSolveFactory.hpp @@ -0,0 +1,511 @@ +// @HEADER +// **************************************************************************** +// Tempus: Copyright (2017) Sandia Corporation +// +// Distributed under BSD 3-clause license (See accompanying file Copyright.txt) +// **************************************************************************** +// @HEADER + +#ifndef Thyra_ReuseLinearOpWithSolveFactory_hpp +#define Thyra_ReuseLinearOpWithSolveFactory_hpp + +#include "Thyra_LinearOpWithSolveFactoryBase.hpp" +#include "Thyra_ReusePreconditionerFactory.hpp" + +namespace Thyra { + +/** \brief A LinearOpWithSolveFactory that is designed to reuse an already + * created/initialized preconditioner. + */ +template +class ReuseLinearOpWithSolveFactory + : virtual public LinearOpWithSolveFactoryBase +{ +public: + + /** @name Overridden from Constructors/Initializers/Accessors */ + //@{ + + /** \brief Construct to uninitialized. */ + ReuseLinearOpWithSolveFactory() {} + + /** \brief Initialize given a single non-const LOWSFB object. + * + * \param lowsf [in,persisting] The LOWSFB object that will be used to + * create the LOWSB object for the diagonal blocks. + * + * Preconditions:
    + *
  • !is_null(lowsf) + *
+ * + */ + void nonconstInitialize( + const RCP > &lowsf, + const RCP > &prec + ); + + + /** \brief Initialize given a single const LOWSFB object. + * + * \param lowsf [in,persisting] The LOWSFB object that will be used to + * create the LOWSB object for the diagonal blocks. + * + * Preconditions:
    + *
  • !is_null(lowsf) + *
+ * + */ + void initialize( + const RCP > &lowsf, + const RCP > &prec + ); + + /** \brief . */ + RCP > getUnderlyingLOWSF(); + + /** \brief . */ + RCP > getUnderlyingLOWSF() const; + + /** \brief . */ + RCP > getUnderlyingPreconditioner(); + + /** \brief . */ + RCP > getUnderlyingPreconditioner() const; + + //@} + + /** \name Overridden from Teuchos::Describable. */ + //@{ + + /** \brief . */ + std::string description() const; + + //@} + + /** @name Overridden from ParameterListAcceptor (simple forwarding functions) */ + //@{ + + /** \brief . */ + void setParameterList(RCP const& paramList); + /** \brief . */ + RCP getNonconstParameterList(); + /** \brief . */ + RCP unsetParameterList(); + /** \brief . */ + RCP getParameterList() const; + /** \brief . */ + RCP getValidParameters() const; + + //@} + + /** \name Overridden from LinearOpWithSolveFactoyBase */ + //@{ + + /** \brief returns false. */ + virtual bool acceptsPreconditionerFactory() const; + + /** \brief Throws exception. */ + virtual void setPreconditionerFactory( + const RCP > &precFactory, + const std::string &precFactoryName + ); + + /** \brief Returns null . */ + virtual RCP > + getPreconditionerFactory() const; + + /** \brief Throws exception. */ + virtual void unsetPreconditionerFactory( + RCP > *precFactory, + std::string *precFactoryName + ); + + /** \brief . */ + virtual bool isCompatible( + const LinearOpSourceBase &fwdOpSrc + ) const; + + /** \brief . */ + virtual RCP > createOp() const; + + /** \brief . */ + virtual void initializeOp( + const RCP > &fwdOpSrc, + LinearOpWithSolveBase *Op, + const ESupportSolveUse supportSolveUse + ) const; + + /** \brief . */ + virtual void initializeAndReuseOp( + const RCP > &fwdOpSrc, + LinearOpWithSolveBase *Op + ) const; + + /** \brief . */ + virtual void uninitializeOp( + LinearOpWithSolveBase *Op, + RCP > *fwdOpSrc, + RCP > *prec, + RCP > *approxFwdOpSrc, + ESupportSolveUse *supportSolveUse + ) const; + + /** \brief . */ + virtual bool supportsPreconditionerInputType( + const EPreconditionerInputType precOpType + ) const; + + /** \brief . */ + virtual void initializePreconditionedOp( + const RCP > &fwdOpSrc, + const RCP > &prec, + LinearOpWithSolveBase *Op, + const ESupportSolveUse supportSolveUse + ) const; + + /** \brief . */ + virtual void initializeApproxPreconditionedOp( + const RCP > &fwdOpSrc, + const RCP > &approxFwdOpSrc, + LinearOpWithSolveBase *Op, + const ESupportSolveUse supportSolveUse + ) const; + + //@} + +protected: + + /** \brief Overridden from Teuchos::VerboseObjectBase */ + //@{ + + /** \brief . */ + void informUpdatedVerbosityState() const; + + //@} + +private: + + typedef Teuchos::ConstNonconstObjectContainer > LOWSF_t; + + LOWSF_t lowsf_; + RCP< PreconditionerBase > prec_; + +}; + +/** \brief Nonmember constructor. + * + * \releates ReuseLinearOpWithSolveFactory + */ +template +RCP > +nonconstReuseLinearOpWithSolveFactory( + const RCP > &lowsf, + const RCP > &prec + ) +{ + RCP > rlowsf = + Teuchos::rcp(new ReuseLinearOpWithSolveFactory); + rlowsf->nonconstInitialize(lowsf, prec); + return rlowsf; +} + +/** \brief Nonmember constructor. + * + * \releates ReuseLinearOpWithSolveFactory + */ +template +RCP > +reuseLinearOpWithSolveFactory( + const RCP > &lowsf, + const RCP > &prec + ) +{ + RCP > rlowsf = + Teuchos::rcp(new ReuseLinearOpWithSolveFactory); + rlowsf->initialize(lowsf, prec); + return rlowsf; +} + +// Overridden from Constructors/Initializers/Accessors + +template +void +ReuseLinearOpWithSolveFactory:: +nonconstInitialize( + const RCP > &lowsf, + const RCP > &prec + ) +{ +#ifdef TEUCHOS_DEBUG + TEUCHOS_TEST_FOR_EXCEPT(is_null(lowsf)); + TEUCHOS_TEST_FOR_EXCEPT(is_null(prec)); +#endif + lowsf_.initialize(lowsf); + prec_ = prec; +} + +template +void +ReuseLinearOpWithSolveFactory:: +initialize( + const RCP > &lowsf, + const RCP > &prec + ) +{ +#ifdef TEUCHOS_DEBUG + TEUCHOS_TEST_FOR_EXCEPT(is_null(lowsf)); + TEUCHOS_TEST_FOR_EXCEPT(is_null(prec)); +#endif + lowsf_.initialize(lowsf); + prec_ = prec; +} + +template +RCP > +ReuseLinearOpWithSolveFactory:: +getUnderlyingLOWSF() +{ + return lowsf_.getNonconstObj(); +} + +template +RCP > +ReuseLinearOpWithSolveFactory:: +getUnderlyingLOWSF() const +{ + return lowsf_.getConstObj(); +} + +template +RCP > +ReuseLinearOpWithSolveFactory:: +getUnderlyingPreconditioner() +{ + return prec_; +} + +template +RCP > +ReuseLinearOpWithSolveFactory:: +getUnderlyingPreconditioner() const +{ + return prec_; +} + +// Overridden from Teuchos::Describable + +template +std::string +ReuseLinearOpWithSolveFactory:: +description() const +{ + std::ostringstream oss; + oss << this->Teuchos::Describable::description() + << "{" + << "lowsf="; + if (!is_null(lowsf_.getConstObj())) + oss << lowsf_.getConstObj()->description(); + else + oss << "NULL"; + oss << std::endl + << "prec="; + if (!is_null(prec_)) + oss << prec_->description(); + else + oss << "NULL"; + oss << "}"; + return oss.str(); +} + +// Overridden from ParameterListAcceptor + +template +void +ReuseLinearOpWithSolveFactory:: +setParameterList( + RCP const& paramList + ) +{ + lowsf_.getNonconstObj()->setParameterList(paramList); +} + +template +RCP +ReuseLinearOpWithSolveFactory:: +getNonconstParameterList() +{ + return lowsf_.getNonconstObj()->getNonconstParameterList(); +} + +template +RCP +ReuseLinearOpWithSolveFactory:: +unsetParameterList() +{ + return lowsf_.getNonconstObj()->unsetParameterList(); +} + +template +RCP +ReuseLinearOpWithSolveFactory:: +getParameterList() const +{ + return lowsf_.getConstObj()->getParameterList(); +} + +template +RCP +ReuseLinearOpWithSolveFactory:: +getValidParameters() const +{ + return lowsf_.getConstObj()->getValidParameters(); +} + +// Overridden from LinearOpWithSolveFactoyBase + +template +bool +ReuseLinearOpWithSolveFactory:: +acceptsPreconditionerFactory() const +{ + return false; +} + +template +void +ReuseLinearOpWithSolveFactory:: +setPreconditionerFactory( + const RCP > &precFactory, + const std::string &precFactoryName + ) +{ +} + +template +RCP > +ReuseLinearOpWithSolveFactory:: +getPreconditionerFactory() const +{ + return Thyra::reusePreconditionerFactory(prec_); +} + +template +void ReuseLinearOpWithSolveFactory:: +unsetPreconditionerFactory( + RCP > *precFactory, + std::string *precFactoryName + ) +{ +} + +template +bool +ReuseLinearOpWithSolveFactory:: +isCompatible( + const LinearOpSourceBase &fwdOpSrc + ) const +{ + return lowsf_.getConstObj()->isCompatible(fwdOpSrc); +} + +template +RCP > +ReuseLinearOpWithSolveFactory:: +createOp() const +{ + return lowsf_.getConstObj()->createOp(); +} + +template +void +ReuseLinearOpWithSolveFactory:: +initializeOp( + const RCP > &fwdOpSrc, + LinearOpWithSolveBase *Op, + const ESupportSolveUse supportSolveUse + ) const +{ + lowsf_.getConstObj()->initializeOp(fwdOpSrc, Op, supportSolveUse); +} + +template +void +ReuseLinearOpWithSolveFactory:: +initializeAndReuseOp( + const RCP > &fwdOpSrc, + LinearOpWithSolveBase *Op + ) const +{ + lowsf_.getConstObj()->initializeAndReuseOp(fwdOpSrc, Op); +} + +template +void +ReuseLinearOpWithSolveFactory:: +uninitializeOp( + LinearOpWithSolveBase *Op, + RCP > *fwdOpSrc, + RCP > *prec, + RCP > *approxFwdOpSrc, + ESupportSolveUse *supportSolveUse + ) const +{ + lowsf_.getConstObj()->uninitializeOp(Op, fwdOpSrc, prec, approxFwdOpSrc, + supportSolveUse); +} + +template +bool +ReuseLinearOpWithSolveFactory:: +supportsPreconditionerInputType( + const EPreconditionerInputType precOpType + ) const +{ + return lowsf_.getConstObj()->supportsPreconditionerInputType(precOpType); +} + +template +void +ReuseLinearOpWithSolveFactory:: +initializePreconditionedOp( + const RCP > &fwdOpSrc, + const RCP > &prec, + LinearOpWithSolveBase *Op, + const ESupportSolveUse supportSolveUse + ) const +{ + lowsf_.getConstObj()->initializePreconditionedOp(fwdOpSrc, prec, Op, + supportSolveUse); +} + +template +void +ReuseLinearOpWithSolveFactory:: +initializeApproxPreconditionedOp( + const RCP > &fwdOpSrc, + const RCP > &approxFwdOpSrc, + LinearOpWithSolveBase *Op, + const ESupportSolveUse supportSolveUse + ) const +{ + lowsf_.getConstObj()->initializeApproxPreconditionedOp(fwdOpSrc, + approxFwdOpSrc, + Op, + supportSolveUse); +} + +// protected + +template +void +ReuseLinearOpWithSolveFactory:: +informUpdatedVerbosityState() const +{ + lowsf_.getConstObj()->setVerbLevel(this->getVerbLevel()); + lowsf_.getConstObj()->setOStream(this->getOStream()); +} + +} // namespace Thyra + + +#endif diff --git a/packages/tempus/src/Thyra_ReusePreconditionerFactory.hpp b/packages/tempus/src/Thyra_ReusePreconditionerFactory.hpp new file mode 100644 index 000000000000..cb12146f94f3 --- /dev/null +++ b/packages/tempus/src/Thyra_ReusePreconditionerFactory.hpp @@ -0,0 +1,179 @@ +// @HEADER +// **************************************************************************** +// Tempus: Copyright (2017) Sandia Corporation +// +// Distributed under BSD 3-clause license (See accompanying file Copyright.txt) +// **************************************************************************** +// @HEADER + +#ifndef Thyra_ReusePreconditionerFactory_hpp +#define Thyra_ReusePreconditionerFactory_hpp + +#include "Thyra_PreconditionerFactoryBase.hpp" + +namespace Thyra { + +/** \brief Concrete PreconditionerFactoryBase subclass that + * just returns an already created/initialized preconditioner object. + */ +template +class ReusePreconditionerFactory + : virtual public PreconditionerFactoryBase +{ +public: + + /** @name Constructors/initializers/accessors */ + //@{ + + /** \brief Construct to uninitialized. */ + ReusePreconditionerFactory() {} + + /** \brief . */ + void initialize( + const RCP > &prec + ) { +#ifdef TEUCHOS_DEBUG + TEUCHOS_TEST_FOR_EXCEPT(is_null(prec)); +#endif + prec_ = prec; + } + + /** \brief . */ + RCP > + getNonconstPreconditioner() { return prec_; } + + /** \brief . */ + RCP > + getPreconditioner() const { return prec_; } + + /** \brief . */ + void uninitialize() { + prec_ = Teuchos::null; + } + + /** \name Overridden from Teuchos::Describable. */ + //@{ + + /** \brief . */ + std::string description() const + { + std::ostringstream oss; + oss << this->Teuchos::Describable::description() + << "{" + << "prec="; + if (!is_null(prec_)) + oss << prec_->description(); + else + oss << "NULL"; + oss << "}"; + return oss.str(); + } + + //@} + + /** @name Overridden from ParameterListAcceptor (simple forwarding functions) */ + //@{ + + /** \brief . */ + void setParameterList(RCP const& paramList) + { + } + + /** \brief . */ + RCP getNonconstParameterList() + { + return Teuchos::null; + } + + /** \brief . */ + RCP unsetParameterList() + { + return Teuchos::null; + } + + /** \brief . */ + RCP getParameterList() const + { + return Teuchos::null; + } + + /** \brief . */ + RCP getValidParameters() const + { + return rcp(new ParameterList); + } + + //@} + + //@} + + /** @name Overridden from PreconditionerFactoryBase */ + //@{ + + /** \brief . */ + bool isCompatible(const LinearOpSourceBase &fwdOpSrc) const + { return false; } + + /** \brief . */ + RCP > createPrec() const + { return prec_; } + + /** \brief . */ + void initializePrec( + const RCP > &fwdOpSrc, + PreconditionerBase *precOp, + const ESupportSolveUse supportSolveUse = SUPPORT_SOLVE_UNSPECIFIED + ) const + { + } + + /** \brief . */ + void uninitializePrec( + PreconditionerBase *precOp, + RCP > *fwdOpSrc = NULL, + ESupportSolveUse *supportSolveUse = NULL + ) const + { + } + + //@} + +private: + + // ////////////////////////////// + // Private data members + + RCP< PreconditionerBase > prec_; + +}; + +/** \brief Nonmember constructor function. + * + * \relates ReusePreconditionerFactory + */ +template +RCP > +reusePreconditionerFactory() +{ + return Teuchos::rcp(new ReusePreconditionerFactory()); +} + +/** \brief Nonmember constructor function. + * + * \relates ReusePreconditionerFactory + */ +template +RCP > +reusePreconditionerFactory( + const RCP > &prec + ) +{ + RCP > + fac = Teuchos::rcp(new ReusePreconditionerFactory()); + fac->initialize(prec); + return fac; +} + +} // end namespace Thyra + +#endif diff --git a/packages/tempus/test/BDF2/CMakeLists.txt b/packages/tempus/test/BDF2/CMakeLists.txt index a4ef9271368c..0052c603e761 100644 --- a/packages/tempus/test/BDF2/CMakeLists.txt +++ b/packages/tempus/test/BDF2/CMakeLists.txt @@ -10,12 +10,12 @@ INCLUDE_DIRECTORIES(REQUIRED_DURING_INSTALLATION_TESTING TRIBITS_ADD_EXECUTABLE( BDF2 - SOURCES Tempus_BDF2Test.cpp ${TEUCHOS_STD_UNIT_TEST_MAIN} + SOURCES Tempus_BDF2Test.cpp ${TEUCHOS_STD_UNIT_TEST_MAIN} TESTONLYLIBS tempus_test_models - COMM serial mpi + COMM serial mpi ) TRIBITS_COPY_FILES_TO_BINARY_DIR(Test_BDF2_CopyFiles - DEST_FILES Tempus_BDF2_SinCos.xml Tempus_BDF2_CDR.xml Tempus_BDF2_VanDerPol.xml + DEST_FILES Tempus_BDF2_SinCos.xml Tempus_BDF2_CDR.xml Tempus_BDF2_VanDerPol.xml Tempus_BDF2_SteadyQuadratic.xml EXEDEPS BDF2 ) diff --git a/packages/tempus/test/BDF2/Tempus_BDF2Test.cpp b/packages/tempus/test/BDF2/Tempus_BDF2Test.cpp index ac06be5d0375..2d1086b82399 100644 --- a/packages/tempus/test/BDF2/Tempus_BDF2Test.cpp +++ b/packages/tempus/test/BDF2/Tempus_BDF2Test.cpp @@ -9,13 +9,17 @@ #include "Teuchos_UnitTestHarness.hpp" #include "Teuchos_XMLParameterListHelpers.hpp" #include "Teuchos_TimeMonitor.hpp" +#include "Teuchos_DefaultComm.hpp" #include "Tempus_config.hpp" #include "Tempus_IntegratorBasic.hpp" +#include "Tempus_IntegratorForwardSensitivity.hpp" +#include "Tempus_IntegratorPseudoTransientForwardSensitivity.hpp" #include "../TestModels/SinCosModel.hpp" #include "../TestModels/CDR_Model.hpp" #include "../TestModels/VanDerPolModel.hpp" +#include "../TestModels/SteadyQuadraticModel.hpp" #include "../TestUtils/Tempus_ConvergenceTestUtils.hpp" #include "Stratimikos_DefaultLinearSolverBuilder.hpp" @@ -191,6 +195,251 @@ TEUCHOS_UNIT_TEST(BDF2, SinCos) ftmp.close(); } +// ************************************************************ +// ************************************************************ +void test_sincos_fsa(const bool use_combined_method, + const bool use_dfdp_as_tangent, + Teuchos::FancyOStream &out, bool &success) +{ + std::vector StepSize; + std::vector ErrorNorm; + const int nTimeStepSizes = 7; + double dt = 0.2; + double order = 0.0; + Teuchos::RCP > comm = + Teuchos::DefaultComm::getComm(); + Teuchos::RCP my_out = + Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout)); + my_out->setProcRankAndSize(comm->getRank(), comm->getSize()); + my_out->setOutputToRootOnly(0); + for (int n=0; n pList = + getParametersFromXmlFile("Tempus_BDF2_SinCos.xml"); + + // Setup the SinCosModel + RCP scm_pl = sublist(pList, "SinCosModel", true); + scm_pl->set("Use DfDp as Tangent", use_dfdp_as_tangent); + RCP > model = + Teuchos::rcp(new SinCosModel(scm_pl)); + + dt /= 2; + + // Setup sensitivities + RCP pl = sublist(pList, "Tempus", true); + ParameterList& sens_pl = pl->sublist("Sensitivities"); + if (use_combined_method) + sens_pl.set("Sensitivity Method", "Combined"); + else { + sens_pl.set("Sensitivity Method", "Staggered"); + sens_pl.set("Reuse State Linear Solver", true); + } + sens_pl.set("Use DfDp as Tangent", use_dfdp_as_tangent); + + // Setup the Integrator and reset initial time step + pl->sublist("Default Integrator") + .sublist("Time Step Control").set("Initial Time Step", dt); + RCP > integrator = + Tempus::integratorForwardSensitivity(pl, model); + order = integrator->getStepper()->getOrder(); + + // Initial Conditions + double t0 = pl->sublist("Default Integrator") + .sublist("Time Step Control").get("Initial Time"); + RCP > x0 = + model->getExactSolution(t0).get_x(); + const int num_param = model->get_p_space(0)->dim(); + RCP > DxDp0 = + Thyra::createMembers(model->get_x_space(), num_param); + for (int i=0; icol(i).ptr(), + *(model->getExactSensSolution(i, t0).get_x())); + integrator->setInitialState(t0, x0, Teuchos::null, Teuchos::null, + DxDp0, Teuchos::null, Teuchos::null); + + // Integrate to timeMax + bool integratorStatus = integrator->advanceTime(); + TEST_ASSERT(integratorStatus) + + // Test if at 'Final Time' + double time = integrator->getTime(); + double timeFinal =pl->sublist("Default Integrator") + .sublist("Time Step Control").get("Final Time"); + TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14); + + // Time-integrated solution and the exact solution + RCP > x = integrator->getX(); + RCP > DxDp = integrator->getDxDp(); + RCP > x_exact = + model->getExactSolution(time).get_x(); + RCP > DxDp_exact = + Thyra::createMembers(model->get_x_space(), num_param); + for (int i=0; icol(i).ptr(), + *(model->getExactSensSolution(i, time).get_x())); + + // Plot sample solution and exact solution + if (comm->getRank() == 0 && n == nTimeStepSizes-1) { + typedef Thyra::DefaultMultiVectorProductVector DMVPV; + + std::ofstream ftmp("Tempus_BDF2_SinCos_Sens.dat"); + RCP > solutionHistory = + integrator->getSolutionHistory(); + RCP< Thyra::MultiVectorBase > DxDp_exact_plot = + Thyra::createMembers(model->get_x_space(), num_param); + for (int i=0; igetNumStates(); i++) { + RCP > solutionState = (*solutionHistory)[i]; + double time = solutionState->getTime(); + RCP x_prod_plot = + Teuchos::rcp_dynamic_cast(solutionState->getX()); + RCP > x_plot = + x_prod_plot->getMultiVector()->col(0); + RCP > DxDp_plot = + x_prod_plot->getMultiVector()->subView(Teuchos::Range1D(1,num_param)); + RCP > x_exact_plot = + model->getExactSolution(time).get_x(); + for (int j=0; jcol(j).ptr(), + *(model->getExactSensSolution(j, time).get_x())); + ftmp << std::fixed << std::setprecision(7) + << time + << std::setw(11) << get_ele(*(x_plot), 0) + << std::setw(11) << get_ele(*(x_plot), 1); + for (int j=0; jcol(j)), 0) + << std::setw(11) << get_ele(*(DxDp_plot->col(j)), 1); + ftmp << std::setw(11) << get_ele(*(x_exact_plot), 0) + << std::setw(11) << get_ele(*(x_exact_plot), 1); + for (int j=0; jcol(j)), 0) + << std::setw(11) << get_ele(*(DxDp_exact_plot->col(j)), 1); + ftmp << std::endl; + } + ftmp.close(); + } + + // Calculate the error + RCP > xdiff = x->clone_v(); + RCP > DxDpdiff = DxDp->clone_mv(); + Thyra::V_StVpStV(xdiff.ptr(), 1.0, *x_exact, -1.0, *(x)); + Thyra::V_VmV(DxDpdiff.ptr(), *DxDp_exact, *DxDp); + StepSize.push_back(dt); + double L2norm = Thyra::norm_2(*xdiff); + L2norm *= L2norm; + Teuchos::Array L2norm_DxDp(num_param); + Thyra::norms_2(*DxDpdiff, L2norm_DxDp()); + for (int i=0; igetRank() == 0) { + std::ofstream ftmp("Tempus_BDF2_SinCos_Sens-Error.dat"); + double error0 = 0.8*ErrorNorm[0]; + for (int n=0; n pList = + getParametersFromXmlFile("Tempus_BDF2_SteadyQuadratic.xml"); + + // Setup the SteadyQuadraticModel + RCP scm_pl = sublist(pList, "SteadyQuadraticModel", true); + scm_pl->set("Use DfDp as Tangent", use_dfdp_as_tangent); + RCP > model = + Teuchos::rcp(new SteadyQuadraticModel(scm_pl)); + + // Setup sensitivities + RCP pl = sublist(pList, "Tempus", true); + ParameterList& sens_pl = pl->sublist("Sensitivities"); + sens_pl.set("Use DfDp as Tangent", use_dfdp_as_tangent); + sens_pl.set("Reuse State Linear Solver", true); + sens_pl.set("Force W Update", true); // Have to do this because for this + // model the solver seems to be overwriting the matrix + + // Setup the Integrator + RCP > integrator = + Tempus::integratorPseudoTransientForwardSensitivity(pl, model); + + // Integrate to timeMax + bool integratorStatus = integrator->advanceTime(); + TEST_ASSERT(integratorStatus); + + // Test if at 'Final Time' + double time = integrator->getTime(); + double timeFinal =pl->sublist("Default Integrator") + .sublist("Time Step Control").get("Final Time"); + TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14); + + // Time-integrated solution and the exact solution + RCP > x_vec = integrator->getX(); + RCP > DxDp_vec = integrator->getDxDp(); + const double x = Thyra::get_ele(*x_vec, 0); + const double dxdb = Thyra::get_ele(*(DxDp_vec->col(0)), 0); + const double x_exact = model->getSteadyStateSolution(); + const double dxdb_exact = model->getSteadyStateSolutionSensitivity(); + + TEST_FLOATING_EQUALITY( x, x_exact, 1.0e-6 ); + TEST_FLOATING_EQUALITY( dxdb, dxdb_exact, 1.0e-6 ); +} + +TEUCHOS_UNIT_TEST(BDF2, SteadyQuadratic_PseudoTransient_FSA) +{ + test_pseudotransient_fsa(false, out, success); +} + +TEUCHOS_UNIT_TEST(BDF2, SteadyQuadratic_PseudoTransient_FSA_Tangent) +{ + test_pseudotransient_fsa(true, out, success); +} + // ************************************************************ // ************************************************************ TEUCHOS_UNIT_TEST(BDF2, CDR) diff --git a/packages/tempus/test/BDF2/Tempus_BDF2_SinCos.xml b/packages/tempus/test/BDF2/Tempus_BDF2_SinCos.xml index 77cb0e37f266..36ece233440a 100644 --- a/packages/tempus/test/BDF2/Tempus_BDF2_SinCos.xml +++ b/packages/tempus/test/BDF2/Tempus_BDF2_SinCos.xml @@ -1,6 +1,6 @@ - + diff --git a/packages/tempus/test/BDF2/Tempus_BDF2_SteadyQuadratic.xml b/packages/tempus/test/BDF2/Tempus_BDF2_SteadyQuadratic.xml new file mode 100644 index 000000000000..f0a46fc9ba11 --- /dev/null +++ b/packages/tempus/test/BDF2/Tempus_BDF2_SteadyQuadratic.xml @@ -0,0 +1,104 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/packages/tempus/test/BackwardEuler/CMakeLists.txt b/packages/tempus/test/BackwardEuler/CMakeLists.txt index 02b272ef3ea4..b9c310e34668 100644 --- a/packages/tempus/test/BackwardEuler/CMakeLists.txt +++ b/packages/tempus/test/BackwardEuler/CMakeLists.txt @@ -9,6 +9,6 @@ TRIBITS_ADD_EXECUTABLE_AND_TEST( ) TRIBITS_COPY_FILES_TO_BINARY_DIR(Test_BackwardEuler_CopyFiles - DEST_FILES Tempus_BackwardEuler_SinCos.xml Tempus_BackwardEuler_CDR.xml Tempus_BackwardEuler_VanDerPol.xml + DEST_FILES Tempus_BackwardEuler_SinCos.xml Tempus_BackwardEuler_CDR.xml Tempus_BackwardEuler_VanDerPol.xml Tempus_BackwardEuler_SteadyQuadratic.xml EXEDEPS BackwardEuler ) diff --git a/packages/tempus/test/BackwardEuler/Tempus_BackwardEulerTest.cpp b/packages/tempus/test/BackwardEuler/Tempus_BackwardEulerTest.cpp index b640993797c2..68b70480f79c 100644 --- a/packages/tempus/test/BackwardEuler/Tempus_BackwardEulerTest.cpp +++ b/packages/tempus/test/BackwardEuler/Tempus_BackwardEulerTest.cpp @@ -9,17 +9,22 @@ #include "Teuchos_UnitTestHarness.hpp" #include "Teuchos_XMLParameterListHelpers.hpp" #include "Teuchos_TimeMonitor.hpp" +#include "Teuchos_DefaultComm.hpp" #include "Tempus_config.hpp" #include "Tempus_IntegratorBasic.hpp" +#include "Tempus_IntegratorForwardSensitivity.hpp" +#include "Tempus_IntegratorPseudoTransientForwardSensitivity.hpp" #include "../TestModels/SinCosModel.hpp" #include "../TestModels/CDR_Model.hpp" #include "../TestModels/VanDerPolModel.hpp" +#include "../TestModels/SteadyQuadraticModel.hpp" #include "../TestUtils/Tempus_ConvergenceTestUtils.hpp" #include "Stratimikos_DefaultLinearSolverBuilder.hpp" #include "Thyra_LinearOpWithSolveFactoryHelpers.hpp" +#include "Thyra_DefaultMultiVectorProductVector.hpp" #ifdef Tempus_ENABLE_MPI #include "Epetra_MpiComm.h" @@ -191,6 +196,251 @@ TEUCHOS_UNIT_TEST(BackwardEuler, SinCos) ftmp.close(); } +// ************************************************************ +// ************************************************************ +void test_sincos_fsa(const bool use_combined_method, + const bool use_dfdp_as_tangent, + Teuchos::FancyOStream &out, bool &success) +{ + std::vector StepSize; + std::vector ErrorNorm; + const int nTimeStepSizes = 7; + double dt = 0.2; + double order = 0.0; + Teuchos::RCP > comm = + Teuchos::DefaultComm::getComm(); + Teuchos::RCP my_out = + Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout)); + my_out->setProcRankAndSize(comm->getRank(), comm->getSize()); + my_out->setOutputToRootOnly(0); + for (int n=0; n pList = + getParametersFromXmlFile("Tempus_BackwardEuler_SinCos.xml"); + + // Setup the SinCosModel + RCP scm_pl = sublist(pList, "SinCosModel", true); + scm_pl->set("Use DfDp as Tangent", use_dfdp_as_tangent); + RCP > model = + Teuchos::rcp(new SinCosModel(scm_pl)); + + dt /= 2; + + // Setup sensitivities + RCP pl = sublist(pList, "Tempus", true); + ParameterList& sens_pl = pl->sublist("Sensitivities"); + if (use_combined_method) + sens_pl.set("Sensitivity Method", "Combined"); + else { + sens_pl.set("Sensitivity Method", "Staggered"); + sens_pl.set("Reuse State Linear Solver", true); + } + sens_pl.set("Use DfDp as Tangent", use_dfdp_as_tangent); + + // Setup the Integrator and reset initial time step + pl->sublist("Default Integrator") + .sublist("Time Step Control").set("Initial Time Step", dt); + RCP > integrator = + Tempus::integratorForwardSensitivity(pl, model); + order = integrator->getStepper()->getOrder(); + + // Initial Conditions + double t0 = pl->sublist("Default Integrator") + .sublist("Time Step Control").get("Initial Time"); + RCP > x0 = + model->getExactSolution(t0).get_x(); + const int num_param = model->get_p_space(0)->dim(); + RCP > DxDp0 = + Thyra::createMembers(model->get_x_space(), num_param); + for (int i=0; icol(i).ptr(), + *(model->getExactSensSolution(i, t0).get_x())); + integrator->setInitialState(t0, x0, Teuchos::null, Teuchos::null, + DxDp0, Teuchos::null, Teuchos::null); + + // Integrate to timeMax + bool integratorStatus = integrator->advanceTime(); + TEST_ASSERT(integratorStatus) + + // Test if at 'Final Time' + double time = integrator->getTime(); + double timeFinal =pl->sublist("Default Integrator") + .sublist("Time Step Control").get("Final Time"); + TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14); + + // Time-integrated solution and the exact solution + RCP > x = integrator->getX(); + RCP > DxDp = integrator->getDxDp(); + RCP > x_exact = + model->getExactSolution(time).get_x(); + RCP > DxDp_exact = + Thyra::createMembers(model->get_x_space(), num_param); + for (int i=0; icol(i).ptr(), + *(model->getExactSensSolution(i, time).get_x())); + + // Plot sample solution and exact solution + if (comm->getRank() == 0 && n == nTimeStepSizes-1) { + typedef Thyra::DefaultMultiVectorProductVector DMVPV; + + std::ofstream ftmp("Tempus_BackwardEuler_SinCos_Sens.dat"); + RCP > solutionHistory = + integrator->getSolutionHistory(); + RCP< Thyra::MultiVectorBase > DxDp_exact_plot = + Thyra::createMembers(model->get_x_space(), num_param); + for (int i=0; igetNumStates(); i++) { + RCP > solutionState = (*solutionHistory)[i]; + double time = solutionState->getTime(); + RCP x_prod_plot = + Teuchos::rcp_dynamic_cast(solutionState->getX()); + RCP > x_plot = + x_prod_plot->getMultiVector()->col(0); + RCP > DxDp_plot = + x_prod_plot->getMultiVector()->subView(Teuchos::Range1D(1,num_param)); + RCP > x_exact_plot = + model->getExactSolution(time).get_x(); + for (int j=0; jcol(j).ptr(), + *(model->getExactSensSolution(j, time).get_x())); + ftmp << std::fixed << std::setprecision(7) + << time + << std::setw(11) << get_ele(*(x_plot), 0) + << std::setw(11) << get_ele(*(x_plot), 1); + for (int j=0; jcol(j)), 0) + << std::setw(11) << get_ele(*(DxDp_plot->col(j)), 1); + ftmp << std::setw(11) << get_ele(*(x_exact_plot), 0) + << std::setw(11) << get_ele(*(x_exact_plot), 1); + for (int j=0; jcol(j)), 0) + << std::setw(11) << get_ele(*(DxDp_exact_plot->col(j)), 1); + ftmp << std::endl; + } + ftmp.close(); + } + + // Calculate the error + RCP > xdiff = x->clone_v(); + RCP > DxDpdiff = DxDp->clone_mv(); + Thyra::V_StVpStV(xdiff.ptr(), 1.0, *x_exact, -1.0, *(x)); + Thyra::V_VmV(DxDpdiff.ptr(), *DxDp_exact, *DxDp); + StepSize.push_back(dt); + double L2norm = Thyra::norm_2(*xdiff); + L2norm *= L2norm; + Teuchos::Array L2norm_DxDp(num_param); + Thyra::norms_2(*DxDpdiff, L2norm_DxDp()); + for (int i=0; igetRank() == 0) { + std::ofstream ftmp("Tempus_BackwardEuler_SinCos_Sens-Error.dat"); + double error0 = 0.8*ErrorNorm[0]; + for (int n=0; n pList = + getParametersFromXmlFile("Tempus_BackwardEuler_SteadyQuadratic.xml"); + + // Setup the SteadyQuadraticModel + RCP scm_pl = sublist(pList, "SteadyQuadraticModel", true); + scm_pl->set("Use DfDp as Tangent", use_dfdp_as_tangent); + RCP > model = + Teuchos::rcp(new SteadyQuadraticModel(scm_pl)); + + // Setup sensitivities + RCP pl = sublist(pList, "Tempus", true); + ParameterList& sens_pl = pl->sublist("Sensitivities"); + sens_pl.set("Use DfDp as Tangent", use_dfdp_as_tangent); + sens_pl.set("Reuse State Linear Solver", true); + sens_pl.set("Force W Update", true); // Have to do this because for this + // model the solver seems to be overwriting the matrix + + // Setup the Integrator + RCP > integrator = + Tempus::integratorPseudoTransientForwardSensitivity(pl, model); + + // Integrate to timeMax + bool integratorStatus = integrator->advanceTime(); + TEST_ASSERT(integratorStatus); + + // Test if at 'Final Time' + double time = integrator->getTime(); + double timeFinal =pl->sublist("Default Integrator") + .sublist("Time Step Control").get("Final Time"); + TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14); + + // Time-integrated solution and the exact solution + RCP > x_vec = integrator->getX(); + RCP > DxDp_vec = integrator->getDxDp(); + const double x = Thyra::get_ele(*x_vec, 0); + const double dxdb = Thyra::get_ele(*(DxDp_vec->col(0)), 0); + const double x_exact = model->getSteadyStateSolution(); + const double dxdb_exact = model->getSteadyStateSolutionSensitivity(); + + TEST_FLOATING_EQUALITY( x, x_exact, 1.0e-6 ); + TEST_FLOATING_EQUALITY( dxdb, dxdb_exact, 1.0e-6 ); +} + +TEUCHOS_UNIT_TEST(BackwardEuler, SteadyQuadratic_PseudoTransient_FSA) +{ + test_pseudotransient_fsa(false, out, success); +} + +TEUCHOS_UNIT_TEST(BackwardEuler, SteadyQuadratic_PseudoTransient_FSA_Tangent) +{ + test_pseudotransient_fsa(true, out, success); +} + // ************************************************************ // ************************************************************ TEUCHOS_UNIT_TEST(BackwardEuler, CDR) diff --git a/packages/tempus/test/BackwardEuler/Tempus_BackwardEuler_SinCos.xml b/packages/tempus/test/BackwardEuler/Tempus_BackwardEuler_SinCos.xml index 30627b6d235f..aadc27948347 100644 --- a/packages/tempus/test/BackwardEuler/Tempus_BackwardEuler_SinCos.xml +++ b/packages/tempus/test/BackwardEuler/Tempus_BackwardEuler_SinCos.xml @@ -1,6 +1,6 @@ - + diff --git a/packages/tempus/test/BackwardEuler/Tempus_BackwardEuler_SteadyQuadratic.xml b/packages/tempus/test/BackwardEuler/Tempus_BackwardEuler_SteadyQuadratic.xml new file mode 100644 index 000000000000..d756e2b3a9c4 --- /dev/null +++ b/packages/tempus/test/BackwardEuler/Tempus_BackwardEuler_SteadyQuadratic.xml @@ -0,0 +1,104 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/packages/tempus/test/ExplicitRK/CMakeLists.txt b/packages/tempus/test/ExplicitRK/CMakeLists.txt index 8751d186bbd3..eba313d32b8b 100644 --- a/packages/tempus/test/ExplicitRK/CMakeLists.txt +++ b/packages/tempus/test/ExplicitRK/CMakeLists.txt @@ -9,6 +9,6 @@ TRIBITS_ADD_EXECUTABLE_AND_TEST( ) TRIBITS_COPY_FILES_TO_BINARY_DIR(Test_ExplicitRK_CopyFiles - DEST_FILES Tempus_ExplicitRK_SinCos.xml + DEST_FILES Tempus_ExplicitRK_SinCos.xml Tempus_ExplicitRK_SteadyQuadratic.xml EXEDEPS ExplicitRK ) diff --git a/packages/tempus/test/ExplicitRK/Tempus_ExplicitRKTest.cpp b/packages/tempus/test/ExplicitRK/Tempus_ExplicitRKTest.cpp index 8c1848a8bb44..8fcaa46155b2 100644 --- a/packages/tempus/test/ExplicitRK/Tempus_ExplicitRKTest.cpp +++ b/packages/tempus/test/ExplicitRK/Tempus_ExplicitRKTest.cpp @@ -9,12 +9,18 @@ #include "Teuchos_UnitTestHarness.hpp" #include "Teuchos_XMLParameterListHelpers.hpp" #include "Teuchos_TimeMonitor.hpp" +#include "Teuchos_DefaultComm.hpp" #include "Thyra_VectorStdOps.hpp" #include "Tempus_IntegratorBasic.hpp" +#include "Tempus_IntegratorForwardSensitivity.hpp" +#include "Tempus_IntegratorPseudoTransientForwardSensitivity.hpp" + +#include "Thyra_DefaultMultiVectorProductVector.hpp" #include "../TestModels/SinCosModel.hpp" +#include "../TestModels/SteadyQuadraticModel.hpp" #include "../TestUtils/Tempus_ConvergenceTestUtils.hpp" #include @@ -251,5 +257,279 @@ TEUCHOS_UNIT_TEST(ExplicitRK, SinCos) Teuchos::TimeMonitor::summarize(); } +// ************************************************************ +// ************************************************************ +void test_sincos_fsa(const bool use_combined_method, + const bool use_dfdp_as_tangent, + Teuchos::FancyOStream &out, bool &success) +{ + std::vector RKMethods; + RKMethods.push_back("RK Forward Euler"); + RKMethods.push_back("RK Explicit 4 Stage"); + RKMethods.push_back("RK Explicit 3/8 Rule"); + RKMethods.push_back("RK Explicit 4 Stage 3rd order by Runge"); + RKMethods.push_back("RK Explicit 5 Stage 3rd order by Kinnmark and Gray"); + RKMethods.push_back("RK Explicit 3 Stage 3rd order"); + RKMethods.push_back("RK Explicit 3 Stage 3rd order TVD"); + RKMethods.push_back("RK Explicit 3 Stage 3rd order by Heun"); + RKMethods.push_back("RK Explicit 2 Stage 2nd order by Runge"); + RKMethods.push_back("RK Explicit Trapezoidal"); + RKMethods.push_back("General ERK"); + std::vector RKMethodErrors; + RKMethodErrors.push_back(0.183799); + RKMethodErrors.push_back(6.88637e-06); + RKMethodErrors.push_back(6.88637e-06); + RKMethodErrors.push_back(0.000264154); + RKMethodErrors.push_back(5.22798e-05); + RKMethodErrors.push_back(0.000261896); + RKMethodErrors.push_back(0.000261896); + RKMethodErrors.push_back(0.000261896); + RKMethodErrors.push_back(0.00934377); + RKMethodErrors.push_back(0.00934377); + RKMethodErrors.push_back(6.88637e-06); + + Teuchos::RCP > comm = + Teuchos::DefaultComm::getComm(); + Teuchos::RCP my_out = + Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout)); + my_out->setProcRankAndSize(comm->getRank(), comm->getSize()); + my_out->setOutputToRootOnly(0); + + for(std::vector::size_type m = 0; m != RKMethods.size(); m++) { + + std::string RKMethod_ = RKMethods[m]; + std::replace(RKMethod_.begin(), RKMethod_.end(), ' ', '_'); + std::replace(RKMethod_.begin(), RKMethod_.end(), '/', '.'); + std::vector StepSize; + std::vector ErrorNorm; + const int nTimeStepSizes = 7; + double dt = 0.2; + double order = 0.0; + for (int n=0; n pList = + getParametersFromXmlFile("Tempus_ExplicitRK_SinCos.xml"); + + // Setup the SinCosModel + RCP scm_pl = sublist(pList, "SinCosModel", true); + scm_pl->set("Use DfDp as Tangent", use_dfdp_as_tangent); + RCP > model = + Teuchos::rcp(new SinCosModel(scm_pl)); + + // Set the Stepper + RCP pl = sublist(pList, "Tempus", true); + if (RKMethods[m] == "General ERK") { + pl->sublist("Demo Integrator").set("Stepper Name", "Demo Stepper 2"); + } else { + pl->sublist("Demo Stepper").set("Stepper Type", RKMethods[m]); + } + + + dt /= 2; + + // Setup sensitivities + ParameterList& sens_pl = pl->sublist("Sensitivities"); + if (use_combined_method) + sens_pl.set("Sensitivity Method", "Combined"); + else + sens_pl.set("Sensitivity Method", "Staggered"); + sens_pl.set("Use DfDp as Tangent", use_dfdp_as_tangent); + + // Setup the Integrator and reset initial time step + pl->sublist("Demo Integrator") + .sublist("Time Step Control").set("Initial Time Step", dt); + RCP > integrator = + Tempus::integratorForwardSensitivity(pl, model); + order = integrator->getStepper()->getOrder(); + + // Initial Conditions + double t0 = pl->sublist("Demo Integrator") + .sublist("Time Step Control").get("Initial Time"); + // RCP > x0 = + // model->getExactSolution(t0).get_x()->clone_v(); + RCP > x0 = + model->getNominalValues().get_x()->clone_v(); + const int num_param = model->get_p_space(0)->dim(); + RCP > DxDp0 = + Thyra::createMembers(model->get_x_space(), num_param); + for (int i=0; icol(i).ptr(), + *(model->getExactSensSolution(i, t0).get_x())); + integrator->setInitialState(t0, x0, Teuchos::null, Teuchos::null, + DxDp0, Teuchos::null, Teuchos::null); + + // Integrate to timeMax + bool integratorStatus = integrator->advanceTime(); + TEST_ASSERT(integratorStatus) + + // Test if at 'Final Time' + double time = integrator->getTime(); + double timeFinal = pl->sublist("Demo Integrator") + .sublist("Time Step Control").get("Final Time"); + TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14); + + // Time-integrated solution and the exact solution + RCP > x = integrator->getX(); + RCP > DxDp = integrator->getDxDp(); + RCP > x_exact = + model->getExactSolution(time).get_x(); + RCP > DxDp_exact = + Thyra::createMembers(model->get_x_space(), num_param); + for (int i=0; icol(i).ptr(), + *(model->getExactSensSolution(i, time).get_x())); + + // Plot sample solution and exact solution + if (comm->getRank() == 0 && n == nTimeStepSizes-1) { + typedef Thyra::DefaultMultiVectorProductVector DMVPV; + + std::ofstream ftmp("Tempus_"+RKMethod_+"_SinCos_Sens.dat"); + RCP > solutionHistory = + integrator->getSolutionHistory(); + RCP< Thyra::MultiVectorBase > DxDp_exact_plot = + Thyra::createMembers(model->get_x_space(), num_param); + for (int i=0; igetNumStates(); i++) { + RCP > solutionState = + (*solutionHistory)[i]; + double time = solutionState->getTime(); + RCP x_prod_plot = + Teuchos::rcp_dynamic_cast(solutionState->getX()); + RCP > x_plot = + x_prod_plot->getMultiVector()->col(0); + RCP > DxDp_plot = + x_prod_plot->getMultiVector()->subView(Teuchos::Range1D(1,num_param)); + RCP > x_exact_plot = + model->getExactSolution(time).get_x(); + for (int j=0; jcol(j).ptr(), + *(model->getExactSensSolution(j, time).get_x())); + ftmp << std::fixed << std::setprecision(7) + << time + << std::setw(11) << get_ele(*(x_plot), 0) + << std::setw(11) << get_ele(*(x_plot), 1); + for (int j=0; jcol(j)), 0) + << std::setw(11) << get_ele(*(DxDp_plot->col(j)), 1); + ftmp << std::setw(11) << get_ele(*(x_exact_plot), 0) + << std::setw(11) << get_ele(*(x_exact_plot), 1); + for (int j=0; jcol(j)), 0) + << std::setw(11) << get_ele(*(DxDp_exact_plot->col(j)), 1); + ftmp << std::endl; + } + ftmp.close(); + } + + // Calculate the error + RCP > xdiff = x->clone_v(); + RCP > DxDpdiff = DxDp->clone_mv(); + Thyra::V_StVpStV(xdiff.ptr(), 1.0, *x_exact, -1.0, *(x)); + Thyra::V_VmV(DxDpdiff.ptr(), *DxDp_exact, *DxDp); + StepSize.push_back(dt); + double L2norm = Thyra::norm_2(*xdiff); + L2norm *= L2norm; + Teuchos::Array L2norm_DxDp(num_param); + Thyra::norms_2(*DxDpdiff, L2norm_DxDp()); + for (int i=0; igetRank() == 0) { + std::ofstream ftmp("Tempus_"+RKMethod_+"_SinCos_Sens-Error.dat"); + double error0 = 0.8*ErrorNorm[0]; + for (int n=0; n pList = + getParametersFromXmlFile("Tempus_ExplicitRK_SteadyQuadratic.xml"); + + // Setup the SteadyQuadraticModel + RCP scm_pl = sublist(pList, "SteadyQuadraticModel", true); + scm_pl->set("Use DfDp as Tangent", use_dfdp_as_tangent); + RCP > model = + Teuchos::rcp(new SteadyQuadraticModel(scm_pl)); + + // Setup sensitivities + RCP pl = sublist(pList, "Tempus", true); + ParameterList& sens_pl = pl->sublist("Sensitivities"); + sens_pl.set("Use DfDp as Tangent", use_dfdp_as_tangent); + + // Setup the Integrator + RCP > integrator = + Tempus::integratorPseudoTransientForwardSensitivity(pl, model); + + // Integrate to timeMax + bool integratorStatus = integrator->advanceTime(); + TEST_ASSERT(integratorStatus); + + // Test if at 'Final Time' + double time = integrator->getTime(); + double timeFinal = pl->sublist("Demo Integrator") + .sublist("Time Step Control").get("Final Time"); + TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14); + + // Time-integrated solution and the exact solution + RCP > x_vec = integrator->getX(); + RCP > DxDp_vec = integrator->getDxDp(); + const double x = Thyra::get_ele(*x_vec, 0); + const double dxdb = Thyra::get_ele(*(DxDp_vec->col(0)), 0); + const double x_exact = model->getSteadyStateSolution(); + const double dxdb_exact = model->getSteadyStateSolutionSensitivity(); + + TEST_FLOATING_EQUALITY( x, x_exact, 1.0e-6 ); + TEST_FLOATING_EQUALITY( dxdb, dxdb_exact, 1.0e-6 ); +} + +TEUCHOS_UNIT_TEST(ExplicitRK, SteadyQuadratic_PseudoTransient_FSA) +{ + test_pseudotransient_fsa(false, out, success); +} + +TEUCHOS_UNIT_TEST(ExplicitRK, SteadyQuadratic_PseudoTransient_FSA_Tangent) +{ + test_pseudotransient_fsa(true, out, success); +} } // namespace Tempus_Test diff --git a/packages/tempus/test/ExplicitRK/Tempus_ExplicitRK_SinCos.xml b/packages/tempus/test/ExplicitRK/Tempus_ExplicitRK_SinCos.xml index 52610a9dca20..cd47a3da799e 100644 --- a/packages/tempus/test/ExplicitRK/Tempus_ExplicitRK_SinCos.xml +++ b/packages/tempus/test/ExplicitRK/Tempus_ExplicitRK_SinCos.xml @@ -1,6 +1,6 @@ - + diff --git a/packages/tempus/test/ExplicitRK/Tempus_ExplicitRK_SteadyQuadratic.xml b/packages/tempus/test/ExplicitRK/Tempus_ExplicitRK_SteadyQuadratic.xml new file mode 100644 index 000000000000..35f8aa74c1ec --- /dev/null +++ b/packages/tempus/test/ExplicitRK/Tempus_ExplicitRK_SteadyQuadratic.xml @@ -0,0 +1,59 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/packages/tempus/test/TestModels/CMakeLists.txt b/packages/tempus/test/TestModels/CMakeLists.txt index dd978eeea206..164b678a3f35 100644 --- a/packages/tempus/test/TestModels/CMakeLists.txt +++ b/packages/tempus/test/TestModels/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}) diff --git a/packages/tempus/test/TestModels/SinCosModel_decl.hpp b/packages/tempus/test/TestModels/SinCosModel_decl.hpp index 25ee0597ac6f..b728f59e9c3c 100644 --- a/packages/tempus/test/TestModels/SinCosModel_decl.hpp +++ b/packages/tempus/test/TestModels/SinCosModel_decl.hpp @@ -149,6 +149,7 @@ class SinCosModel int ng_; ///< Number of elements in this observation function (1) bool haveIC_; ///< false => no nominal values are provided (default=true) bool acceptModelParams_; ///< Changes inArgs to require parameters + bool useDfDpAsTangent_; ///< Treat DfDp OutArg as tangent (df/dx*dx/dp+df/dp) mutable bool isInitialized_; mutable Thyra::ModelEvaluatorBase::InArgs inArgs_; mutable Thyra::ModelEvaluatorBase::OutArgs outArgs_; @@ -157,6 +158,7 @@ class SinCosModel Teuchos::RCP > f_space_; Teuchos::RCP > p_space_; Teuchos::RCP > g_space_; + Teuchos::RCP > DxDp_space_; // Parameters for the model: x_0(t) = a + b*sin(f*t+phi) // x_1(t) = b*f*cos(f*t+phi) diff --git a/packages/tempus/test/TestModels/SinCosModel_impl.hpp b/packages/tempus/test/TestModels/SinCosModel_impl.hpp index 1c7abc3654f3..11cef18b5ee6 100644 --- a/packages/tempus/test/TestModels/SinCosModel_impl.hpp +++ b/packages/tempus/test/TestModels/SinCosModel_impl.hpp @@ -18,6 +18,7 @@ #include "Thyra_DefaultMultiVectorLinearOpWithSolve.hpp" #include "Thyra_DefaultLinearOpSource.hpp" #include "Thyra_VectorStdOps.hpp" +#include "Thyra_DefaultMultiVectorProductVector.hpp" #include @@ -30,11 +31,12 @@ SinCosModel(Teuchos::RCP pList_) { isInitialized_ = false; dim_ = 2; - Np_ = 1; // Number of parameter vectors (1) + Np_ = 3; // Number of parameter vectors (p, dx/dp, dx_dot/dp) np_ = 3; // Number of parameters in this vector (3) Ng_ = 1; // Number of observation functions (1) ng_ = 1; // Number of elements in this observation function (1) acceptModelParams_ = false; + useDfDpAsTangent_ = false; haveIC_ = true; a_ = 0.0; f_ = 1.0; @@ -51,6 +53,9 @@ SinCosModel(Teuchos::RCP pList_) g_space_ = Thyra::defaultSpmdVectorSpace(ng_); setParameterList(pList_); + + // Create DxDp product space + DxDp_space_ = Thyra::multiVectorProductVectorSpace(x_space_, np_); } template @@ -275,8 +280,11 @@ evalModelImpl( const Thyra::ModelEvaluatorBase::OutArgs &outArgs ) const { + typedef Thyra::DefaultMultiVectorProductVector DMVPV; using Teuchos::RCP; using Thyra::VectorBase; + using Thyra::MultiVectorBase; + using Teuchos::rcp_dynamic_cast; TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error, "Error, setupInOutArgs_ must be called first!\n"); @@ -296,6 +304,16 @@ evalModelImpl( L = p_in_view[2]; } + RCP > DxDp_in, DxdotDp_in; + if (acceptModelParams_) { + if (inArgs.get_p(1) != Teuchos::null) + DxDp_in = + rcp_dynamic_cast(inArgs.get_p(1))->getMultiVector(); + if (inArgs.get_p(2) != Teuchos::null) + DxdotDp_in = + rcp_dynamic_cast(inArgs.get_p(2))->getMultiVector(); + } + Scalar beta = inArgs.get_beta(); const RCP > f_out = outArgs.get_f(); @@ -332,6 +350,17 @@ evalModelImpl( DfDp_out_view(1,0) = (f/L)*(f/L); DfDp_out_view(1,1) = (2.0*f/(L*L))*(a-x_in_view[0]); DfDp_out_view(1,2) = -(2.0*f*f/(L*L*L))*(a-x_in_view[0]); + + // Compute df/dp + (df/dx) * (dx/dp) + if (useDfDpAsTangent_ && !is_null(DxDp_in)) { + Thyra::ConstDetachedMultiVectorView DxDp( *DxDp_in ); + DfDp_out_view(0,0) += DxDp(1,0); + DfDp_out_view(0,1) += DxDp(1,1); + DfDp_out_view(0,2) += DxDp(1,2); + DfDp_out_view(1,0) += -(f/L)*(f/L) * DxDp(0,0); + DfDp_out_view(1,1) += -(f/L)*(f/L) * DxDp(0,1); + DfDp_out_view(1,2) += -(f/L)*(f/L) * DxDp(0,2); + } } } else { @@ -363,6 +392,26 @@ evalModelImpl( DfDp_out_view(1,0) = -(f/L)*(f/L); DfDp_out_view(1,1) = -(2.0*f/(L*L))*(a-x_in_view[0]); DfDp_out_view(1,2) = +(2.0*f*f/(L*L*L))*(a-x_in_view[0]); + + // Compute df/dp + (df/dx_dot) * (dx_dot/dp) + (df/dx) * (dx/dp) + if (useDfDpAsTangent_ && !is_null(DxdotDp_in)) { + Thyra::ConstDetachedMultiVectorView DxdotDp( *DxdotDp_in ); + DfDp_out_view(0,0) += DxdotDp(0,0); + DfDp_out_view(0,1) += DxdotDp(0,1); + DfDp_out_view(0,2) += DxdotDp(0,2); + DfDp_out_view(1,0) += DxdotDp(1,0); + DfDp_out_view(1,1) += DxdotDp(1,1); + DfDp_out_view(1,2) += DxdotDp(1,2); + } + if (useDfDpAsTangent_ && !is_null(DxDp_in)) { + Thyra::ConstDetachedMultiVectorView DxDp( *DxDp_in ); + DfDp_out_view(0,0) += -DxDp(1,0); + DfDp_out_view(0,1) += -DxDp(1,1); + DfDp_out_view(0,2) += -DxDp(1,2); + DfDp_out_view(1,0) += (f/L)*(f/L) * DxDp(0,0); + DfDp_out_view(1,1) += (f/L)*(f/L) * DxDp(0,1); + DfDp_out_view(1,2) += (f/L)*(f/L) * DxDp(0,2); + } } } } @@ -376,7 +425,11 @@ get_p_space(int l) const return Teuchos::null; } TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE( l, 0, Np_ ); - return p_space_; + if (l == 0) + return p_space_; + else if (l == 1 || l == 2) + return DxDp_space_; + return Teuchos::null; } template @@ -390,9 +443,15 @@ get_p_names(int l) const TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE( l, 0, Np_ ); Teuchos::RCP > p_strings = Teuchos::rcp(new Teuchos::Array()); - p_strings->push_back("Model Coefficient: a"); - p_strings->push_back("Model Coefficient: f"); - p_strings->push_back("Model Coefficient: L"); + if (l == 0) { + p_strings->push_back("Model Coefficient: a"); + p_strings->push_back("Model Coefficient: f"); + p_strings->push_back("Model Coefficient: L"); + } + else if (l == 1) + p_strings->push_back("DxDp"); + else if (l == 2) + p_strings->push_back("Dx_dotDp"); return p_strings; } @@ -494,6 +553,7 @@ setParameterList(Teuchos::RCP const& paramList) Teuchos::RCP pl = this->getMyNonconstParamList(); bool acceptModelParams = get(*pl,"Accept model parameters"); bool haveIC = get(*pl,"Provide nominal values"); + bool useDfDpAsTangent = get(*pl, "Use DfDp as Tangent"); if ( (acceptModelParams != acceptModelParams_) || (haveIC != haveIC_) ) { @@ -501,6 +561,7 @@ setParameterList(Teuchos::RCP const& paramList) } acceptModelParams_ = acceptModelParams; haveIC_ = haveIC; + useDfDpAsTangent_ = useDfDpAsTangent; a_ = get(*pl,"Coeff a"); f_ = get(*pl,"Coeff f"); L_ = get(*pl,"Coeff L"); @@ -521,6 +582,7 @@ getValidParameters() const Teuchos::RCP pl = Teuchos::parameterList(); pl->set("Accept model parameters", false); pl->set("Provide nominal values", true); + pl->set("Use DfDp as Tangent", false); Teuchos::setDoubleParameter( "Coeff a", 0.0, "Coefficient a in model", &*pl); Teuchos::setDoubleParameter( diff --git a/packages/tempus/test/TestModels/SteadyQuadraticModel.cpp b/packages/tempus/test/TestModels/SteadyQuadraticModel.cpp new file mode 100644 index 000000000000..9def3e70b7af --- /dev/null +++ b/packages/tempus/test/TestModels/SteadyQuadraticModel.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 "SteadyQuadraticModel.hpp" +#include "SteadyQuadraticModel_impl.hpp" + +namespace Tempus_Test { + TEMPUS_INSTANTIATE_TEMPLATE_CLASS(SteadyQuadraticModel) +} + +#endif diff --git a/packages/tempus/test/TestModels/SteadyQuadraticModel_decl.hpp b/packages/tempus/test/TestModels/SteadyQuadraticModel_decl.hpp new file mode 100644 index 000000000000..95ce058673a7 --- /dev/null +++ b/packages/tempus/test/TestModels/SteadyQuadraticModel_decl.hpp @@ -0,0 +1,105 @@ +// @HEADER +// **************************************************************************** +// Tempus: Copyright (2017) Sandia Corporation +// +// Distributed under BSD 3-clause license (See accompanying file Copyright.txt) +// **************************************************************************** +// @HEADER + +#ifndef TEMPUS_TEST_STEADY_QUADRATIC_MODEL_DECL_HPP +#define TEMPUS_TEST_STEADY_QUADRATIC_MODEL_DECL_HPP + +#include "Thyra_ModelEvaluator.hpp" // Interface +#include "Thyra_StateFuncModelEvaluatorBase.hpp" // Implementation + +#include "Teuchos_ParameterListAcceptorDefaultBase.hpp" +#include "Teuchos_ParameterList.hpp" + +namespace Tempus_Test { + +/** \brief Simple quadratic equation with a stable steady-state. + * This is a simple differential equation + * \f[ + * \mathbf{\dot{x}}=\mathbf{x}^2 - b^2 + * \f] + * which has steady state solutions \f$\mathbf{x} = \pm b\f$. The solution + * \f$\mathbf{x} = b\f$ is stable if \f$b < 0\f$ and the solution + * \f$\mathbf{x} = -b\f$ is stable if \f$b > 0\f$. This model is used to + * test pseudo-transient sensitivity analysis methods. + */ +template +class SteadyQuadraticModel + : public Thyra::StateFuncModelEvaluatorBase, + public Teuchos::ParameterListAcceptorDefaultBase +{ + public: + + // Constructor + SteadyQuadraticModel(Teuchos::RCP pList = Teuchos::null); + + // Exact solution + Scalar getSteadyStateSolution() const; + + // Exact sensitivity solution + Scalar getSteadyStateSolutionSensitivity() const; + + /** \name Public functions overridden from ModelEvaluator. */ + //@{ + + Teuchos::RCP > get_x_space() const; + Teuchos::RCP > get_f_space() const; + Thyra::ModelEvaluatorBase::InArgs getNominalValues() const; + Teuchos::RCP > create_W() const; + Teuchos::RCP > create_W_op() const; + Teuchos::RCP > get_W_factory() const; + Thyra::ModelEvaluatorBase::InArgs createInArgs() const; + + Teuchos::RCP > get_p_space(int l) const; + Teuchos::RCP > get_p_names(int l) const; + Teuchos::RCP > get_g_space(int j) const; + + //@} + + /** \name Public functions overridden from ParameterListAcceptor. */ + //@{ + void setParameterList(Teuchos::RCP const& paramList); + Teuchos::RCP getValidParameters() const; + //@} + +private: + + void setupInOutArgs_() const; + + /** \name Private functions overridden from ModelEvaluatorDefaultBase. */ + //@{ + Thyra::ModelEvaluatorBase::OutArgs createOutArgsImpl() const; + void evalModelImpl( + const Thyra::ModelEvaluatorBase::InArgs &inArgs_bar, + const Thyra::ModelEvaluatorBase::OutArgs &outArgs_bar + ) const; + //@} + +private: + int dim_; ///< Number of state unknowns (2) + int Np_; ///< Number of parameter vectors (1) + int np_; ///< Number of parameters in this vector (2) + int Ng_; ///< Number of observation functions (1) + int ng_; ///< Number of elements in this observation function (1) + bool useDfDpAsTangent_; ///< Treat DfDp OutArg as tangent (df/dx*dx/dp+df/dp) + mutable bool isInitialized_; + mutable Thyra::ModelEvaluatorBase::InArgs inArgs_; + mutable Thyra::ModelEvaluatorBase::OutArgs outArgs_; + mutable Thyra::ModelEvaluatorBase::InArgs nominalValues_; + Teuchos::RCP > x_space_; + Teuchos::RCP > f_space_; + Teuchos::RCP > p_space_; + Teuchos::RCP > g_space_; + Teuchos::RCP > DxDp_space_; + + // Parameters for the model + Scalar b_; ///< Model parameter +}; + + +} // namespace Tempus_Test +#endif // TEMPUS_TEST_STEADY_QUADRATIC_MODEL_DECL_HPP diff --git a/packages/tempus/test/TestModels/SteadyQuadraticModel_impl.hpp b/packages/tempus/test/TestModels/SteadyQuadraticModel_impl.hpp new file mode 100644 index 000000000000..82e80d33f13b --- /dev/null +++ b/packages/tempus/test/TestModels/SteadyQuadraticModel_impl.hpp @@ -0,0 +1,420 @@ +// @HEADER +// **************************************************************************** +// Tempus: Copyright (2017) Sandia Corporation +// +// Distributed under BSD 3-clause license (See accompanying file Copyright.txt) +// **************************************************************************** +// @HEADER + +#ifndef TEMPUS_TEST_STEADY_QUADRATIC_MODEL_IMPL_HPP +#define TEMPUS_TEST_STEADY_QUADRATIC_MODEL_IMPL_HPP + +#include "Teuchos_StandardParameterEntryValidators.hpp" + +#include "Thyra_DefaultSpmdVectorSpace.hpp" +#include "Thyra_DetachedVectorView.hpp" +#include "Thyra_DetachedMultiVectorView.hpp" +#include "Thyra_DefaultSerialDenseLinearOpWithSolveFactory.hpp" +#include "Thyra_DefaultMultiVectorLinearOpWithSolve.hpp" +#include "Thyra_DefaultLinearOpSource.hpp" +#include "Thyra_VectorStdOps.hpp" +#include "Thyra_DefaultMultiVectorProductVector.hpp" + +#include + + +namespace Tempus_Test { + +template +SteadyQuadraticModel:: +SteadyQuadraticModel(Teuchos::RCP pList_) +{ + isInitialized_ = false; + dim_ = 1; + Np_ = 3; // Number of parameter vectors (p, dx/dp, dx_dot/dp) + np_ = 1; // Number of parameters in this vector (3) + Ng_ = 1; // Number of observation functions (1) + ng_ = 1; // Number of elements in this observation function (1) + useDfDpAsTangent_ = false; + b_ = 1.0; + + // Create x_space and f_space + x_space_ = Thyra::defaultSpmdVectorSpace(dim_); + f_space_ = Thyra::defaultSpmdVectorSpace(dim_); + // Create p_space and g_space + p_space_ = Thyra::defaultSpmdVectorSpace(np_); + g_space_ = Thyra::defaultSpmdVectorSpace(ng_); + + setParameterList(pList_); + + // Create DxDp product space + DxDp_space_ = Thyra::multiVectorProductVectorSpace(x_space_, np_); +} + +template +Scalar +SteadyQuadraticModel:: +getSteadyStateSolution() const +{ + if (b_ < 0.0) + return b_; + return -b_; +} + +template +Scalar +SteadyQuadraticModel:: +getSteadyStateSolutionSensitivity() const +{ + if (b_ < 0.0) + return 1.0; + return -1.0; +} + +template +Teuchos::RCP > +SteadyQuadraticModel:: +get_x_space() const +{ + return x_space_; +} + + +template +Teuchos::RCP > +SteadyQuadraticModel:: +get_f_space() const +{ + return f_space_; +} + + +template +Thyra::ModelEvaluatorBase::InArgs +SteadyQuadraticModel:: +getNominalValues() const +{ + TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error, + "Error, setupInOutArgs_ must be called first!\n"); + return nominalValues_; +} + + +template +Teuchos::RCP > +SteadyQuadraticModel:: +create_W() const +{ + using Teuchos::RCP; + RCP > W_factory = this->get_W_factory(); + RCP > matrix = this->create_W_op(); + { + // 01/20/09 tscoffe: This is a total hack to provide a full rank matrix to + // linearOpWithSolve because it ends up factoring the matrix during + // initialization, which it really shouldn't do, or I'm doing something + // wrong here. The net effect is that I get exceptions thrown in + // optimized mode due to the matrix being rank deficient unless I do this. + RCP > multivec = + Teuchos::rcp_dynamic_cast >(matrix,true); + { + RCP > vec = Thyra::createMember(x_space_); + { + Thyra::DetachedVectorView vec_view( *vec ); + vec_view[0] = 1.0; + } + V_V(multivec->col(0).ptr(),*vec); + } + } + RCP > W = + Thyra::linearOpWithSolve( + *W_factory, + matrix + ); + return W; +} + +template +Teuchos::RCP > +SteadyQuadraticModel:: +create_W_op() const +{ + Teuchos::RCP > matrix = + Thyra::createMembers(x_space_, dim_); + return(matrix); +} + + +template +Teuchos::RCP > +SteadyQuadraticModel:: +get_W_factory() const +{ + Teuchos::RCP > W_factory = + Thyra::defaultSerialDenseLinearOpWithSolveFactory(); + return W_factory; +} + + +template +Thyra::ModelEvaluatorBase::InArgs +SteadyQuadraticModel:: +createInArgs() const +{ + setupInOutArgs_(); + return inArgs_; +} + + +// Private functions overridden from ModelEvaluatorDefaultBase + + +template +Thyra::ModelEvaluatorBase::OutArgs +SteadyQuadraticModel:: +createOutArgsImpl() const +{ + setupInOutArgs_(); + return outArgs_; +} + + +template +void +SteadyQuadraticModel:: +evalModelImpl( + const Thyra::ModelEvaluatorBase::InArgs &inArgs, + const Thyra::ModelEvaluatorBase::OutArgs &outArgs + ) const +{ + typedef Thyra::DefaultMultiVectorProductVector DMVPV; + using Teuchos::RCP; + using Thyra::VectorBase; + using Thyra::MultiVectorBase; + using Teuchos::rcp_dynamic_cast; + TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error, + "Error, setupInOutArgs_ must be called first!\n"); + + const RCP > x_in = inArgs.get_x().assert_not_null(); + Thyra::ConstDetachedVectorView x_in_view( *x_in ); + + Scalar b = b_; + const RCP > p_in = inArgs.get_p(0); + if (p_in != Teuchos::null) { + Thyra::ConstDetachedVectorView p_in_view( *p_in ); + b = p_in_view[0]; + } + + RCP > DxDp_in, DxdotDp_in; + if (inArgs.get_p(1) != Teuchos::null) + DxDp_in = + rcp_dynamic_cast(inArgs.get_p(1))->getMultiVector(); + if (inArgs.get_p(2) != Teuchos::null) + DxdotDp_in = + rcp_dynamic_cast(inArgs.get_p(2))->getMultiVector(); + + Scalar beta = inArgs.get_beta(); + + const RCP > f_out = outArgs.get_f(); + const RCP > W_out = outArgs.get_W_op(); + RCP > DfDp_out; + Thyra::ModelEvaluatorBase::Derivative DfDp = outArgs.get_DfDp(0); + DfDp_out = DfDp.getMultiVector(); + + if (inArgs.get_x_dot().is_null()) { + + // Evaluate the Explicit ODE f(x,t) [= 0] + if (!is_null(f_out)) { + Thyra::DetachedVectorView f_out_view( *f_out ); + f_out_view[0] = x_in_view[0]*x_in_view[0] - b*b; + } + if (!is_null(W_out)) { + RCP > matrix = + Teuchos::rcp_dynamic_cast >(W_out,true); + Thyra::DetachedMultiVectorView matrix_view( *matrix ); + matrix_view(0,0) = beta*2.0*x_in_view[0]; + } + if (!is_null(DfDp_out)) { + Thyra::DetachedMultiVectorView DfDp_out_view( *DfDp_out ); + DfDp_out_view(0,0) = -2.0*b; + + // Compute df/dp + (df/dx) * (dx/dp) + if (useDfDpAsTangent_ && !is_null(DxDp_in)) { + Thyra::ConstDetachedMultiVectorView DxDp( *DxDp_in ); + DfDp_out_view(0,0) += 2.0*x_in_view[0]*DxDp(0,0); + } + } + } else { + + // Evaluate the implicit ODE f(xdot, x, t) [= 0] + RCP > x_dot_in; + x_dot_in = inArgs.get_x_dot().assert_not_null(); + Scalar alpha = inArgs.get_alpha(); + if (!is_null(f_out)) { + Thyra::DetachedVectorView f_out_view( *f_out ); + Thyra::ConstDetachedVectorView x_dot_in_view( *x_dot_in ); + f_out_view[0] = x_dot_in_view[0] - (x_in_view[0]*x_in_view[0] - b*b); + } + if (!is_null(W_out)) { + RCP > matrix = + Teuchos::rcp_dynamic_cast >(W_out,true); + Thyra::DetachedMultiVectorView matrix_view( *matrix ); + matrix_view(0,0) = alpha - beta*2.0*x_in_view[0]; + } + if (!is_null(DfDp_out)) { + Thyra::DetachedMultiVectorView DfDp_out_view( *DfDp_out ); + DfDp_out_view(0,0) = 2.0*b; + + // Compute df/dp + (df/dx_dot) * (dx_dot/dp) + (df/dx) * (dx/dp) + if (useDfDpAsTangent_ && !is_null(DxdotDp_in)) { + Thyra::ConstDetachedMultiVectorView DxdotDp( *DxdotDp_in ); + DfDp_out_view(0,0) += DxdotDp(0,0); + } + if (useDfDpAsTangent_ && !is_null(DxDp_in)) { + Thyra::ConstDetachedMultiVectorView DxDp( *DxDp_in ); + DfDp_out_view(0,0) += -2.0*x_in_view[0]*DxDp(0,0); + } + } + } +} + +template +Teuchos::RCP > +SteadyQuadraticModel:: +get_p_space(int l) const +{ + TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE( l, 0, Np_ ); + if (l == 0) + return p_space_; + else if (l == 1 || l == 2) + return DxDp_space_; + return Teuchos::null; +} + +template +Teuchos::RCP > +SteadyQuadraticModel:: +get_p_names(int l) const +{ + TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE( l, 0, Np_ ); + Teuchos::RCP > p_strings = + Teuchos::rcp(new Teuchos::Array()); + if (l == 0) { + p_strings->push_back("Model Coefficient: b"); + } + else if (l == 1) + p_strings->push_back("DxDp"); + else if (l == 2) + p_strings->push_back("Dx_dotDp"); + return p_strings; +} + +template +Teuchos::RCP > +SteadyQuadraticModel:: +get_g_space(int j) const +{ + TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE( j, 0, Ng_ ); + return g_space_; +} + +// private + +template +void +SteadyQuadraticModel:: +setupInOutArgs_() const +{ + if (isInitialized_) { + return; + } + + { + // Set up prototypical InArgs + Thyra::ModelEvaluatorBase::InArgsSetup inArgs; + inArgs.setModelEvalDescription(this->description()); + inArgs.setSupports( Thyra::ModelEvaluatorBase::IN_ARG_t ); + inArgs.setSupports( Thyra::ModelEvaluatorBase::IN_ARG_x ); + inArgs.setSupports( Thyra::ModelEvaluatorBase::IN_ARG_beta ); + inArgs.setSupports( Thyra::ModelEvaluatorBase::IN_ARG_x_dot ); + inArgs.setSupports( Thyra::ModelEvaluatorBase::IN_ARG_alpha ); + inArgs.set_Np(Np_); + inArgs_ = inArgs; + } + + { + // Set up prototypical OutArgs + Thyra::ModelEvaluatorBase::OutArgsSetup outArgs; + outArgs.setModelEvalDescription(this->description()); + outArgs.setSupports( Thyra::ModelEvaluatorBase::OUT_ARG_f ); + outArgs.setSupports( Thyra::ModelEvaluatorBase::OUT_ARG_W_op ); + outArgs.set_Np_Ng(Np_,Ng_); + outArgs.setSupports( Thyra::ModelEvaluatorBase::OUT_ARG_DfDp,0, + Thyra::ModelEvaluatorBase::DERIV_MV_BY_COL ); + outArgs_ = outArgs; + } + + // Set up nominal values + nominalValues_ = inArgs_; + nominalValues_.set_t(0.0); + const Teuchos::RCP > x_ic = createMember(x_space_); + { // scope to delete DetachedVectorView + Thyra::DetachedVectorView x_ic_view( *x_ic ); + x_ic_view[0] = 0.0; + } + nominalValues_.set_x(x_ic); + const Teuchos::RCP > p_ic = createMember(p_space_); + { + Thyra::DetachedVectorView p_ic_view( *p_ic ); + p_ic_view[0] = b_; + } + nominalValues_.set_p(0,p_ic); + + const Teuchos::RCP > x_dot_ic = + createMember(x_space_); + { // scope to delete DetachedVectorView + Thyra::DetachedVectorView x_dot_ic_view( *x_dot_ic ); + x_dot_ic_view[0] = 0.0; + } + nominalValues_.set_x_dot(x_dot_ic); + + isInitialized_ = true; + +} + +template +void +SteadyQuadraticModel:: +setParameterList(Teuchos::RCP const& paramList) +{ + using Teuchos::get; + using Teuchos::ParameterList; + Teuchos::RCP tmpPL = + Teuchos::rcp(new ParameterList("SteadyQuadraticModel")); + if (paramList != Teuchos::null) tmpPL = paramList; + tmpPL->validateParametersAndSetDefaults(*this->getValidParameters()); + this->setMyParamList(tmpPL); + Teuchos::RCP pl = this->getMyNonconstParamList(); + useDfDpAsTangent_ = get(*pl, "Use DfDp as Tangent"); + b_ = get(*pl,"Coeff b"); + isInitialized_ = false; // For setup of new in/out args + setupInOutArgs_(); +} + +template +Teuchos::RCP +SteadyQuadraticModel:: +getValidParameters() const +{ + static Teuchos::RCP validPL; + if (is_null(validPL)) { + Teuchos::RCP pl = Teuchos::parameterList(); + pl->set("Use DfDp as Tangent", false); + Teuchos::setDoubleParameter( + "Coeff b", 1.0, "Coefficient b in model", &*pl); + validPL = pl; + } + return validPL; +} + +} // namespace Tempus_Test +#endif // TEMPUS_TEST_STEADY_QUADRATIC_MODEL_IMPL_HPP