From e8e6d6768c2b52d881b236342462322fe396c1a9 Mon Sep 17 00:00:00 2001 From: Eric Phipps Date: Thu, 9 Nov 2017 14:07:19 -0700 Subject: [PATCH] Tempus: Add combined, staggered, and pseudotransient FSA methods This commit adds integrators supporting the combined, staggered, and pseudotransient forward sensitivity analysis methods where the sensitivity equations are solved alongside the forward equations. The different methods refer to how the solves between states/sensitivities are partitioned: * combined: solved together (best for explicit) * staggered: solve state equations and then sensitivity equations for each time step (best for implicit) * pseudotransient: do full state equation integration followed by full sensitivity equation integration (for when time integrator is used as a steady-state solver) In all cases there is an option, "Use DfDp as Tangent", that controls whether the DfDp out-arg is interpreted as the tangent df/dp + df/dx*dx/dp + df/dxdot*dxdot/dp, or not (which is not regular model evaluator behavior). If the model supports it, it is a much simpler and more efficient way to assemble the sensitivity equations. It also supports options to reuse the linear solver from the state equations when solving the sensitivity equations. This saves Jacobian assembly time and preconditioner generation time. Supporting it required one small change to the NOX-Thyra interface to not throw when accessing the preconditioner RCP if it is null (because there is no way to know beforehand if it is null or not until you access it). Tests were added for Backward Euler, Explicit RK, and BDF2 time integration methods for all of the above methods and options. --- packages/nox/src-thyra/NOX_Thyra_Group.C | 2 - packages/tempus/src/CMakeLists.txt | 3 +- ...mbinedForwardSensitivityModelEvaluator.cpp | 19 + ...dForwardSensitivityModelEvaluator_decl.hpp | 136 ++++ ...dForwardSensitivityModelEvaluator_impl.hpp | 413 +++++++++++ .../src/Tempus_IntegratorBasic_decl.hpp | 9 +- .../src/Tempus_IntegratorBasic_impl.hpp | 19 +- .../Tempus_IntegratorForwardSensitivity.cpp | 37 + ...mpus_IntegratorForwardSensitivity_decl.hpp | 259 +++++++ ...mpus_IntegratorForwardSensitivity_impl.hpp | 359 ++++++++++ ...ratorPseudoTransientForwardSensitivity.cpp | 37 + ...PseudoTransientForwardSensitivity_decl.hpp | 190 +++++ ...PseudoTransientForwardSensitivity_impl.hpp | 581 +++++++++++++++ .../src/Tempus_SolutionHistory_decl.hpp | 4 + .../src/Tempus_SolutionStateMetaData_decl.hpp | 2 +- .../src/Tempus_SolutionStateMetaData_impl.hpp | 2 +- .../tempus/src/Tempus_SolutionState_decl.hpp | 10 +- ...ggeredForwardSensitivityModelEvaluator.cpp | 19 + ...dForwardSensitivityModelEvaluator_decl.hpp | 168 +++++ ...dForwardSensitivityModelEvaluator_impl.hpp | 381 ++++++++++ packages/tempus/src/Tempus_Stepper.hpp | 2 + .../tempus/src/Tempus_StepperBDF2_decl.hpp | 2 + .../src/Tempus_StepperBackwardEuler_decl.hpp | 2 + .../tempus/src/Tempus_StepperDIRK_decl.hpp | 2 + .../src/Tempus_StepperExplicitRK_decl.hpp | 2 + .../src/Tempus_StepperForwardEuler_decl.hpp | 2 + .../src/Tempus_StepperHHTAlpha_decl.hpp | 2 + .../Tempus_StepperIMEX_RK_Partition_decl.hpp | 2 + .../tempus/src/Tempus_StepperIMEX_RK_decl.hpp | 2 + .../src/Tempus_StepperLeapfrog_decl.hpp | 2 + ...empus_StepperNewmarkExplicitAForm_decl.hpp | 2 + ...empus_StepperNewmarkImplicitAForm_decl.hpp | 2 + ...empus_StepperNewmarkImplicitDForm_decl.hpp | 2 + ...pus_StepperStaggeredForwardSensitivity.cpp | 19 + ...tepperStaggeredForwardSensitivity_decl.hpp | 143 ++++ ...tepperStaggeredForwardSensitivity_impl.hpp | 425 +++++++++++ .../tempus/src/Thyra_MultiVectorLinearOp.hpp | 359 ++++++++++ ...ra_MultiVectorLinearOpWithSolveFactory.hpp | 668 ++++++++++++++++++ .../src/Thyra_MultiVectorPreconditioner.hpp | 213 ++++++ ...Thyra_MultiVectorPreconditionerFactory.hpp | 274 +++++++ .../Thyra_ReuseLinearOpWithSolveFactory.hpp | 511 ++++++++++++++ .../src/Thyra_ReusePreconditionerFactory.hpp | 179 +++++ packages/tempus/test/BDF2/CMakeLists.txt | 6 +- packages/tempus/test/BDF2/Tempus_BDF2Test.cpp | 249 +++++++ .../tempus/test/BDF2/Tempus_BDF2_SinCos.xml | 2 +- .../test/BDF2/Tempus_BDF2_SteadyQuadratic.xml | 104 +++ .../tempus/test/BackwardEuler/CMakeLists.txt | 2 +- .../Tempus_BackwardEulerTest.cpp | 250 +++++++ .../Tempus_BackwardEuler_SinCos.xml | 2 +- .../Tempus_BackwardEuler_SteadyQuadratic.xml | 104 +++ .../tempus/test/ExplicitRK/CMakeLists.txt | 2 +- .../test/ExplicitRK/Tempus_ExplicitRKTest.cpp | 280 ++++++++ .../ExplicitRK/Tempus_ExplicitRK_SinCos.xml | 2 +- .../Tempus_ExplicitRK_SteadyQuadratic.xml | 59 ++ .../tempus/test/TestModels/CMakeLists.txt | 2 +- .../test/TestModels/SinCosModel_decl.hpp | 2 + .../test/TestModels/SinCosModel_impl.hpp | 72 +- .../test/TestModels/SteadyQuadraticModel.cpp | 19 + .../TestModels/SteadyQuadraticModel_decl.hpp | 105 +++ .../TestModels/SteadyQuadraticModel_impl.hpp | 420 +++++++++++ 60 files changed, 7122 insertions(+), 27 deletions(-) create mode 100644 packages/tempus/src/Tempus_CombinedForwardSensitivityModelEvaluator.cpp create mode 100644 packages/tempus/src/Tempus_CombinedForwardSensitivityModelEvaluator_decl.hpp create mode 100644 packages/tempus/src/Tempus_CombinedForwardSensitivityModelEvaluator_impl.hpp create mode 100644 packages/tempus/src/Tempus_IntegratorForwardSensitivity.cpp create mode 100644 packages/tempus/src/Tempus_IntegratorForwardSensitivity_decl.hpp create mode 100644 packages/tempus/src/Tempus_IntegratorForwardSensitivity_impl.hpp create mode 100644 packages/tempus/src/Tempus_IntegratorPseudoTransientForwardSensitivity.cpp create mode 100644 packages/tempus/src/Tempus_IntegratorPseudoTransientForwardSensitivity_decl.hpp create mode 100644 packages/tempus/src/Tempus_IntegratorPseudoTransientForwardSensitivity_impl.hpp create mode 100644 packages/tempus/src/Tempus_StaggeredForwardSensitivityModelEvaluator.cpp create mode 100644 packages/tempus/src/Tempus_StaggeredForwardSensitivityModelEvaluator_decl.hpp create mode 100644 packages/tempus/src/Tempus_StaggeredForwardSensitivityModelEvaluator_impl.hpp create mode 100644 packages/tempus/src/Tempus_StepperStaggeredForwardSensitivity.cpp create mode 100644 packages/tempus/src/Tempus_StepperStaggeredForwardSensitivity_decl.hpp create mode 100644 packages/tempus/src/Tempus_StepperStaggeredForwardSensitivity_impl.hpp create mode 100644 packages/tempus/src/Thyra_MultiVectorLinearOp.hpp create mode 100644 packages/tempus/src/Thyra_MultiVectorLinearOpWithSolveFactory.hpp create mode 100644 packages/tempus/src/Thyra_MultiVectorPreconditioner.hpp create mode 100644 packages/tempus/src/Thyra_MultiVectorPreconditionerFactory.hpp create mode 100644 packages/tempus/src/Thyra_ReuseLinearOpWithSolveFactory.hpp create mode 100644 packages/tempus/src/Thyra_ReusePreconditionerFactory.hpp create mode 100644 packages/tempus/test/BDF2/Tempus_BDF2_SteadyQuadratic.xml create mode 100644 packages/tempus/test/BackwardEuler/Tempus_BackwardEuler_SteadyQuadratic.xml create mode 100644 packages/tempus/test/ExplicitRK/Tempus_ExplicitRK_SteadyQuadratic.xml create mode 100644 packages/tempus/test/TestModels/SteadyQuadraticModel.cpp create mode 100644 packages/tempus/test/TestModels/SteadyQuadraticModel_decl.hpp create mode 100644 packages/tempus/test/TestModels/SteadyQuadraticModel_impl.hpp 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