9 #ifndef TEMPUS_TEST_VANDERPOL_IMEX_EXPLICITMODEL_IMPL_HPP 
   10 #define TEMPUS_TEST_VANDERPOL_IMEX_EXPLICITMODEL_IMPL_HPP 
   12 #include "Teuchos_StandardParameterEntryValidators.hpp" 
   14 #include "Thyra_DefaultSpmdVectorSpace.hpp" 
   15 #include "Thyra_DetachedVectorView.hpp" 
   16 #include "Thyra_DetachedMultiVectorView.hpp" 
   17 #include "Thyra_DefaultSerialDenseLinearOpWithSolveFactory.hpp" 
   18 #include "Thyra_DefaultMultiVectorLinearOpWithSolve.hpp" 
   19 #include "Thyra_DefaultLinearOpSource.hpp" 
   20 #include "Thyra_VectorStdOps.hpp" 
   21 #include "Thyra_MultiVectorStdOps.hpp" 
   22 #include "Thyra_DefaultMultiVectorProductVector.hpp" 
   23 #include "Thyra_DefaultBlockedLinearOp.hpp" 
   28 namespace Tempus_Test {
 
   30 template<
class Scalar>
 
   33   Teuchos::RCP<Teuchos::ParameterList> pList_, 
bool useProductVector)
 
   35   useProductVector_ = useProductVector;
 
   36   isInitialized_ = 
false;
 
   42   acceptModelParams_ = 
false;
 
   49   if (useProductVector_ == 
false) {
 
   50     x_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(dim_);
 
   51     f_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(dim_);
 
   55     Teuchos::Array<RCP<const Thyra::VectorSpaceBase<Scalar> > > yxSpaceArray;
 
   56     xSpace_ = Thyra::defaultSpmdVectorSpace<Scalar>(1);
 
   57     for (
int i=0; i < dim_; ++i) yxSpaceArray.push_back(xSpace_);
 
   58     x_space_ = Thyra::productVectorSpace<Scalar>(yxSpaceArray);
 
   59     f_space_ = Thyra::productVectorSpace<Scalar>(yxSpaceArray);
 
   63   p_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(np_);
 
   64   g_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(ng_);
 
   65   dxdp_space_ = Thyra::multiVectorProductVectorSpace(x_space_, np_);
 
   67   setParameterList(pList_);
 
   70 template<
class Scalar>
 
   71 Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
 
   79 template<
class Scalar>
 
   80 Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
 
   88 template<
class Scalar>
 
   89 Thyra::ModelEvaluatorBase::InArgs<Scalar>
 
   93   TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
 
   94       "Error, setupInOutArgs_ must be called first!\n");
 
   95   return nominalValues_;
 
   99 template<
class Scalar>
 
  100 Teuchos::RCP<Thyra::LinearOpBase<Scalar> >
 
  104   Teuchos::RCP<Thyra::LinearOpBase<Scalar> > W_op;
 
  105   if (useProductVector_ == 
false)
 
  106     W_op = Thyra::createMembers(x_space_, dim_);
 
  110     Teuchos::RCP<Thyra::LinearOpBase<Scalar> > A00 =
 
  111       Thyra::createMembers(xSpace_, 1);
 
  112     Teuchos::RCP<Thyra::LinearOpBase<Scalar> > A01 =
 
  113       Thyra::createMembers(xSpace_, 1);
 
  114     Teuchos::RCP<Thyra::LinearOpBase<Scalar> > A10 =
 
  115       Thyra::createMembers(xSpace_, 1);
 
  116     Teuchos::RCP<Thyra::LinearOpBase<Scalar> > A11 =
 
  117       Thyra::createMembers(xSpace_, 1);
 
  118     W_op = Thyra::nonconstBlock2x2(A00, A01, A10, A11);
 
  124 template<
class Scalar>
 
  125 Thyra::ModelEvaluatorBase::InArgs<Scalar>
 
  134 template<
class Scalar>
 
  135 Thyra::ModelEvaluatorBase::OutArgs<Scalar>
 
  144 template<
class Scalar>
 
  149   if (isInitialized_) {
 
  155     Thyra::ModelEvaluatorBase::InArgsSetup<Scalar> inArgs;
 
  156     inArgs.setModelEvalDescription(this->description());
 
  157     inArgs.setSupports( Thyra::ModelEvaluatorBase::IN_ARG_x_dot );
 
  158     inArgs.setSupports( Thyra::ModelEvaluatorBase::IN_ARG_x );
 
  159     inArgs.setSupports( Thyra::ModelEvaluatorBase::IN_ARG_t );
 
  160     inArgs.setSupports( Thyra::ModelEvaluatorBase::IN_ARG_alpha );
 
  161     inArgs.setSupports( Thyra::ModelEvaluatorBase::IN_ARG_beta );
 
  162     if (acceptModelParams_) {
 
  170     Thyra::ModelEvaluatorBase::OutArgsSetup<Scalar> outArgs;
 
  171     outArgs.setModelEvalDescription(this->description());
 
  172     outArgs.setSupports( Thyra::ModelEvaluatorBase::OUT_ARG_f );
 
  173     outArgs.setSupports( Thyra::ModelEvaluatorBase::OUT_ARG_W_op );
 
  174     if (acceptModelParams_) {
 
  175       outArgs.set_Np_Ng(Np_,Ng_);
 
  176       outArgs.setSupports( Thyra::ModelEvaluatorBase::OUT_ARG_DfDp,0,
 
  177                            Thyra::ModelEvaluatorBase::DERIV_MV_JACOBIAN_FORM );
 
  183   nominalValues_ = inArgs_;
 
  188     nominalValues_.set_t(t0_ic_);
 
  189     const RCP<Thyra::VectorBase<Scalar> > x_ic = createMember(x_space_);
 
  191       Thyra::DetachedVectorView<Scalar> x_ic_view( *x_ic );
 
  192       x_ic_view[0] = x0_ic_;
 
  193       x_ic_view[1] = x1_ic_;
 
  195     nominalValues_.set_x(x_ic);
 
  197     if (acceptModelParams_) {
 
  198       const RCP<Thyra::VectorBase<Scalar> > p_ic = createMember(p_space_);
 
  200         Thyra::DetachedVectorView<Scalar> p_ic_view( *p_ic );
 
  201         p_ic_view[0] = epsilon_;
 
  203       nominalValues_.set_p(0,p_ic);
 
  205     const RCP<Thyra::VectorBase<Scalar> > x_dot_ic = createMember(x_space_);
 
  207       Thyra::DetachedVectorView<Scalar> x_dot_ic_view( *x_dot_ic );
 
  208       x_dot_ic_view[0] = x1_ic_;
 
  209       x_dot_ic_view[1] = ((1.0-x0_ic_*x0_ic_)*x1_ic_-x0_ic_)/epsilon_;
 
  211     nominalValues_.set_x_dot(x_dot_ic);
 
  214   isInitialized_ = 
true;
 
  217 template<
class Scalar>
 
  223   using Teuchos::ParameterList;
 
  224   Teuchos::RCP<ParameterList> tmpPL = Teuchos::rcp(
new ParameterList(
"VanDerPol_IMEX_ExplicitModel"));
 
  225   if (paramList != Teuchos::null) tmpPL = paramList;
 
  226   tmpPL->validateParametersAndSetDefaults(*this->getValidParameters());
 
  227   this->setMyParamList(tmpPL);
 
  228   Teuchos::RCP<ParameterList> pl = this->getMyNonconstParamList();
 
  229   bool acceptModelParams = get<bool>(*pl,
"Accept model parameters");
 
  230   bool haveIC = get<bool>(*pl,
"Provide nominal values");
 
  231   bool useDfDpAsTangent = get<bool>(*pl, 
"Use DfDp as Tangent");
 
  232   if ( (acceptModelParams != acceptModelParams_) ||
 
  235     isInitialized_ = 
false;
 
  237   acceptModelParams_ = acceptModelParams;
 
  239   useDfDpAsTangent_ = useDfDpAsTangent;
 
  240   epsilon_ = get<Scalar>(*pl,
"Coeff epsilon");
 
  241   x0_ic_ = get<Scalar>(*pl,
"IC x0");
 
  242   x1_ic_ = get<Scalar>(*pl,
"IC x1");
 
  243   t0_ic_ = get<Scalar>(*pl,
"IC t0");
 
  247 template<
class Scalar>
 
  248 Teuchos::RCP<const Teuchos::ParameterList>
 
  252   static Teuchos::RCP<const Teuchos::ParameterList> validPL;
 
  253   if (is_null(validPL)) {
 
  254     Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
 
  255     pl->set(
"Accept model parameters", 
false);
 
  256     pl->set(
"Provide nominal values", 
true);
 
  257     pl->set(
"Use DfDp as Tangent", 
false);
 
  258     Teuchos::setDoubleParameter(
 
  259         "Coeff epsilon", 1.0e-06, 
"Coefficient a in model", &*pl);
 
  260     Teuchos::setDoubleParameter(
 
  261         "IC x0", 2.0, 
"Initial Condition for x0", &*pl);
 
  262     Teuchos::setDoubleParameter(
 
  263         "IC x1", 0.0, 
"Initial Condition for x1", &*pl);
 
  264     Teuchos::setDoubleParameter(
 
  265         "IC t0", 0.0, 
"Initial time t0", &*pl);
 
  271 template<
class Scalar>
 
  272 Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
 
  276   if (!acceptModelParams_) {
 
  277     return Teuchos::null;
 
  279   TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE( l, 0, Np_ );
 
  282   else if (l == 1 || l == 2)
 
  284   return Teuchos::null;
 
  287 template<
class Scalar>
 
  288 Teuchos::RCP<const Teuchos::Array<std::string> >
 
  292   if (!acceptModelParams_) {
 
  293     return Teuchos::null;
 
  295   TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE( l, 0, Np_ );
 
  296   Teuchos::RCP<Teuchos::Array<std::string> > p_strings =
 
  297     Teuchos::rcp(
new Teuchos::Array<std::string>());
 
  299     p_strings->push_back(
"Model Coefficient:  epsilon");
 
  301     p_strings->push_back(
"DxDp");
 
  303     p_strings->push_back(
"Dx_dotDp");
 
  307 template<
class Scalar>
 
  308 Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
 
  312   TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE( j, 0, Ng_ );
 
  316 template<
class Scalar>
 
  320               const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs)
 const 
  322   typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;
 
  324   using Teuchos::rcp_dynamic_cast;
 
  325   TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
 
  326       "Error, setupInOutArgs_ must be called first!\n");
 
  328   const RCP<const Thyra::VectorBase<Scalar> > x_in =
 
  329     inArgs.get_x().assert_not_null();
 
  330   Thyra::ConstDetachedVectorView<Scalar> x_in_view( *x_in );
 
  333   Scalar beta  = inArgs.get_beta();
 
  334   Scalar eps = epsilon_;
 
  335   RCP<const Thyra::MultiVectorBase<Scalar> > dxdp_in, dxdotdp_in;
 
  336   if (acceptModelParams_) {
 
  337     const RCP<const Thyra::VectorBase<Scalar> > p_in =
 
  339     if (p_in != Teuchos::null) {
 
  340       Thyra::ConstDetachedVectorView<Scalar> p_in_view( *p_in );
 
  343     if (inArgs.get_p(1) != Teuchos::null) {
 
  345         rcp_dynamic_cast<
const DMVPV>(inArgs.get_p(1),
true)->getMultiVector();
 
  347     if (inArgs.get_p(2) != Teuchos::null) {
 
  349         rcp_dynamic_cast<
const DMVPV>(inArgs.get_p(2),
true)->getMultiVector();
 
  353   const RCP<Thyra::VectorBase<Scalar> > f_out = outArgs.get_f();
 
  354   const RCP<Thyra::LinearOpBase<Scalar> > W_out = outArgs.get_W_op();
 
  355   RCP<Thyra::MultiVectorBase<Scalar> > DfDp_out;
 
  356   if (acceptModelParams_) {
 
  357     Thyra::ModelEvaluatorBase::Derivative<Scalar> DfDp = outArgs.get_DfDp(0);
 
  358     DfDp_out = DfDp.getMultiVector();
 
  361   if (inArgs.get_x_dot().is_null()) {
 
  364     if (!is_null(f_out)) {
 
  365       Thyra::DetachedVectorView<Scalar> f_out_view( *f_out );
 
  366       f_out_view[0] = x_in_view[1];
 
  367       f_out_view[1] = -x_in_view[0]/eps;
 
  369     if (!is_null(DfDp_out)) {
 
  370       Thyra::DetachedMultiVectorView<Scalar> DfDp_out_view( *DfDp_out );
 
  371       DfDp_out_view(0,0) = 0.0;
 
  372       DfDp_out_view(1,0) = x_in_view[0]/(eps*eps);
 
  375       if (useDfDpAsTangent_ && !is_null(dxdp_in)) {
 
  376         Thyra::ConstDetachedMultiVectorView<Scalar> dxdp( *dxdp_in );
 
  377         DfDp_out_view(0,0) += dxdp(1,0);
 
  378         DfDp_out_view(1,0) += -dxdp(0,0)/eps;
 
  381     if (!is_null(W_out)) {
 
  382       if (useProductVector_ == 
false) {
 
  383         RCP<Thyra::MultiVectorBase<Scalar> > W =
 
  384           Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(W_out,
true);
 
  385         Thyra::DetachedMultiVectorView<Scalar> W_view( *W );
 
  388         W_view(1,0) = -beta/eps;                               
 
  393         RCP<Thyra::DefaultBlockedLinearOp<Scalar> > W_block =
 
  394           Teuchos::rcp_dynamic_cast<Thyra::DefaultBlockedLinearOp<Scalar> >(W_out, 
true);
 
  395         RCP<Thyra::MultiVectorBase<Scalar> > W00 =
 
  396           Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(W_block->getNonconstBlock(0,0),
true);
 
  397         RCP<Thyra::MultiVectorBase<Scalar> > W01 =
 
  398           Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(W_block->getNonconstBlock(0,1),
true);
 
  399         RCP<Thyra::MultiVectorBase<Scalar> > W10 =
 
  400           Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(W_block->getNonconstBlock(1,0),
true);
 
  401         RCP<Thyra::MultiVectorBase<Scalar> > W11 =
 
  402           Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(W_block->getNonconstBlock(1,1),
true);
 
  403         Thyra::DetachedMultiVectorView<Scalar> W00_view( *W00 );
 
  404         Thyra::DetachedMultiVectorView<Scalar> W01_view( *W01 );
 
  405         Thyra::DetachedMultiVectorView<Scalar> W10_view( *W10 );
 
  406         Thyra::DetachedMultiVectorView<Scalar> W11_view( *W11 );
 
  408         W01_view(0,0) = beta;                               
 
  409         W10_view(0,0) = -beta/eps;                          
 
  417     RCP<const Thyra::VectorBase<Scalar> > x_dot_in;
 
  418     x_dot_in = inArgs.get_x_dot().assert_not_null();
 
  419     Scalar alpha = inArgs.get_alpha();
 
  421     if (!is_null(f_out)) {
 
  422       Thyra::DetachedVectorView<Scalar> f_out_view( *f_out );
 
  423       Thyra::ConstDetachedVectorView<Scalar> x_dot_in_view( *x_dot_in );
 
  424       f_out_view[0] = x_dot_in_view[0] - x_in_view[1];
 
  425       f_out_view[1] = x_dot_in_view[1] + x_in_view[0]/eps;
 
  427     if (!is_null(DfDp_out)) {
 
  428       Thyra::DetachedMultiVectorView<Scalar> DfDp_out_view( *DfDp_out );
 
  429       DfDp_out_view(0,0) = 0.0;
 
  430       DfDp_out_view(1,0) = -x_in_view[0]/(eps*eps);
 
  433       if (useDfDpAsTangent_ && !is_null(dxdotdp_in)) {
 
  434         Thyra::ConstDetachedMultiVectorView<Scalar> dxdotdp( *dxdotdp_in );
 
  435         DfDp_out_view(0,0) += dxdotdp(0,0);
 
  436         DfDp_out_view(1,0) += dxdotdp(1,0);
 
  438       if (useDfDpAsTangent_ && !is_null(dxdp_in)) {
 
  439         Thyra::ConstDetachedMultiVectorView<Scalar> dxdp( *dxdp_in );
 
  440         DfDp_out_view(0,0) += -dxdp(1,0);
 
  441         DfDp_out_view(1,0) += dxdp(0,0)/eps;
 
  444     if (!is_null(W_out)) {
 
  445       if (useProductVector_ == 
false) {
 
  446         RCP<Thyra::MultiVectorBase<Scalar> > W =
 
  447           Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(W_out,
true);
 
  448         Thyra::DetachedMultiVectorView<Scalar> W_view( *W );
 
  451         W_view(1,0) = beta/eps;                                
 
  456         RCP<Thyra::DefaultBlockedLinearOp<Scalar> > W_block =
 
  457           Teuchos::rcp_dynamic_cast<Thyra::DefaultBlockedLinearOp<Scalar> >(W_out, 
true);
 
  458         RCP<Thyra::MultiVectorBase<Scalar> > W00 =
 
  459           Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(W_block->getNonconstBlock(0,0),
true);
 
  460         RCP<Thyra::MultiVectorBase<Scalar> > W01 =
 
  461           Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(W_block->getNonconstBlock(0,1),
true);
 
  462         RCP<Thyra::MultiVectorBase<Scalar> > W10 =
 
  463           Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(W_block->getNonconstBlock(1,0),
true);
 
  464         RCP<Thyra::MultiVectorBase<Scalar> > W11 =
 
  465           Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(W_block->getNonconstBlock(1,1),
true);
 
  466         Thyra::DetachedMultiVectorView<Scalar> W00_view( *W00 );
 
  467         Thyra::DetachedMultiVectorView<Scalar> W01_view( *W01 );
 
  468         Thyra::DetachedMultiVectorView<Scalar> W10_view( *W10 );
 
  469         Thyra::DetachedMultiVectorView<Scalar> W11_view( *W11 );
 
  470         W00_view(0,0) = alpha;                               
 
  471         W01_view(0,0) = -beta;                               
 
  472         W10_view(0,0) = beta/eps;                            
 
  473         W11_view(0,0) = alpha;                               
 
  482 #endif // TEMPUS_TEST_VANDERPOL_IMEX_EXPLICITMODEL_IMPL_HPP 
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const 
 
void setupInOutArgs_() const 
 
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_f_space() const 
 
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_x_space() const 
 
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_p_space(int l) const 
 
void setParameterList(Teuchos::RCP< Teuchos::ParameterList > const ¶mList)
 
Teuchos::RCP< Thyra::LinearOpBase< Scalar > > create_W_op() const 
 
Thyra::ModelEvaluatorBase::InArgs< Scalar > getNominalValues() const 
 
Thyra::ModelEvaluatorBase::InArgs< Scalar > createInArgs() const 
 
Thyra::ModelEvaluatorBase::OutArgs< Scalar > createOutArgsImpl() const 
 
VanDerPol_IMEX_ExplicitModel(Teuchos::RCP< Teuchos::ParameterList > pList=Teuchos::null, bool useProductVector=false)
 
void evalModelImpl(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs_bar, const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &outArgs_bar) const 
 
Teuchos::RCP< const Teuchos::Array< std::string > > get_p_names(int l) const 
 
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_g_space(int j) const