69 using Teuchos::rcp_dynamic_cast;
71 this->useImplicitModel_ =
true;
72 this->wrapperImplicitInArgs_ = this->createInArgs();
73 this->wrapperImplicitOutArgs_ = this->createOutArgs();
74 this->useImplicitModel_ =
false;
76 RCP<const Thyra::VectorBase<Scalar> > z =
77 this->explicitModel_->getNominalValues().get_x();
80 TEUCHOS_TEST_FOR_EXCEPTION( !(getIMEXVector(z)->space()->isCompatible(
81 *(this->implicitModel_->get_x_space()))),
83 "Error - WrapperModelEvaluatorPairIMEX_CombinedFSA::initialize()\n"
84 " Explicit and Implicit vector x spaces are incompatible!\n"
85 " Explicit vector x space = " << *(getIMEXVector(z)->space()) <<
"\n"
86 " Implicit vector x space = " << *(this->implicitModel_->get_x_space()) <<
90 const RCP<const DMVPV> z_dmvpv = rcp_dynamic_cast<const DMVPV>(z,
true);
91 const RCP<const Thyra::MultiVectorBase<Scalar> > z_mv =
92 z_dmvpv->getMultiVector();
93 RCP<const PMVB> zPVector = rcp_dynamic_cast<const PMVB>(z_mv);
94 TEUCHOS_TEST_FOR_EXCEPTION( zPVector == Teuchos::null, std::logic_error,
95 "Error - WrapperModelEvaluatorPairPartIMEX_CombinedFSA::initialize()\n"
96 " was given a VectorBase that could not be cast to a\n"
97 " ProductVectorBase!\n");
99 int numBlocks = zPVector->productSpace()->numBlocks();
101 TEUCHOS_TEST_FOR_EXCEPTION( !(0 <= this->numExplicitOnlyBlocks_ &&
102 this->numExplicitOnlyBlocks_ < numBlocks),
104 "Error - WrapperModelEvaluatorPairPartIMEX_CombinedFSA::initialize()\n"
105 "Invalid number of explicit-only blocks = " <<
106 this->numExplicitOnlyBlocks_ <<
"\n"
107 "Should be in the interval [0, numBlocks) = [0, " << numBlocks <<
")\n");
330evalModelImpl(
const Thyra::ModelEvaluatorBase::InArgs<Scalar> & inArgs,
331 const Thyra::ModelEvaluatorBase::OutArgs<Scalar> & outArgs)
const
333 typedef Thyra::ModelEvaluatorBase MEB;
335 using Teuchos::rcp_dynamic_cast;
336 using Teuchos::Range1D;
338 const int p_index = this->parameterIndex_;
343 RCP<const Thyra::VectorBase<Scalar> > x = inArgs.get_x();
344 RCP<Thyra::VectorBase<Scalar> > x_dot =
345 Thyra::createMember(fsaImplicitModel_->get_x_space());
346 this->timeDer_->compute(x, x_dot);
348 MEB::InArgs<Scalar> fsaImplicitInArgs (this->wrapperImplicitInArgs_);
349 MEB::OutArgs<Scalar> fsaImplicitOutArgs(this->wrapperImplicitOutArgs_);
350 fsaImplicitInArgs.set_x(x);
351 fsaImplicitInArgs.set_x_dot(x_dot);
352 for (
int i=0; i<fsaImplicitModel_->Np(); ++i) {
354 if ((inArgs.get_p(i) != Teuchos::null) && (i != p_index))
355 fsaImplicitInArgs.set_p(i, inArgs.get_p(i));
362 RCP<const Thyra::VectorBase<Scalar> > y;
363 RCP<const Thyra::MultiVectorBase<Scalar> > dydp;
364 if (fsaImplicitInArgs.get_p(p_index) != Teuchos::null) {
365 RCP<const Thyra::VectorBase<Scalar> > p =
366 fsaImplicitInArgs.get_p(p_index);
367 RCP<const Thyra::MultiVectorBase<Scalar> > p_mv =
368 rcp_dynamic_cast<const DMVPV>(p,
true)->getMultiVector();
369 const int num_param = p_mv->domain()->dim()-1;
371 dydp = p_mv->subView(Range1D(1,num_param));
372 fsaImplicitInArgs.set_p(p_index, y);
374 if (use_dfdp_as_tangent_) {
375 RCP< const Thyra::VectorBase<Scalar> > dydp_vec =
376 Thyra::multiVectorProductVector(explicit_dydp_prod_space_, dydp);
377 fsaImplicitInArgs.set_p(y_tangent_index_, dydp_vec);
380 fsaImplicitOutArgs.set_f(outArgs.get_f());
381 fsaImplicitOutArgs.set_W_op(outArgs.get_W_op());
383 fsaImplicitModel_->evalModel(fsaImplicitInArgs,fsaImplicitOutArgs);
387 if (!use_dfdp_as_tangent_ && outArgs.get_f() != Teuchos::null) {
388 MEB::InArgs<Scalar> appImplicitInArgs =
389 appImplicitModel_->getNominalValues();
390 RCP< const Thyra::VectorBase<Scalar> > app_x =
391 rcp_dynamic_cast<const DMVPV>(x,
true)->getMultiVector()->col(0);
392 RCP< const Thyra::VectorBase<Scalar> > app_x_dot =
393 rcp_dynamic_cast<const DMVPV>(x_dot,
true)->getMultiVector()->col(0);
394 appImplicitInArgs.set_x(app_x);
395 appImplicitInArgs.set_x_dot(app_x_dot);
396 for (
int i=0; i<appImplicitModel_->Np(); ++i) {
398 appImplicitInArgs.set_p(i, inArgs.get_p(i));
400 appImplicitInArgs.set_p(p_index, y);
401 if (appImplicitInArgs.supports(MEB::IN_ARG_t))
402 appImplicitInArgs.set_t(inArgs.get_t());
403 MEB::OutArgs<Scalar> appImplicitOutArgs =
404 appImplicitModel_->createOutArgs();
405 MEB::DerivativeSupport dfdp_support =
406 appImplicitOutArgs.supports(MEB::OUT_ARG_DfDp, p_index);
407 Thyra::EOpTransp trans = Thyra::NOTRANS;
408 if (dfdp_support.supports(MEB::DERIV_LINEAR_OP)) {
409 if (my_dfdp_op_ == Teuchos::null)
410 my_dfdp_op_ = appImplicitModel_->create_DfDp_op(p_index);
411 appImplicitOutArgs.set_DfDp(p_index,
412 MEB::Derivative<Scalar>(my_dfdp_op_));
413 trans = Thyra::NOTRANS;
415 else if (dfdp_support.supports(MEB::DERIV_MV_JACOBIAN_FORM)) {
416 if (my_dfdp_mv_ == Teuchos::null)
417 my_dfdp_mv_ = Thyra::createMembers(
418 appImplicitModel_->get_f_space(),
419 appImplicitModel_->get_p_space(p_index)->dim());
420 appImplicitOutArgs.set_DfDp(
421 p_index, MEB::Derivative<Scalar>(my_dfdp_mv_,
422 MEB::DERIV_MV_JACOBIAN_FORM));
423 my_dfdp_op_ = my_dfdp_mv_;
424 trans = Thyra::NOTRANS;
426 else if (dfdp_support.supports(MEB::DERIV_MV_GRADIENT_FORM)) {
427 if (my_dfdp_mv_ == Teuchos::null)
428 my_dfdp_mv_ = Thyra::createMembers(
429 appImplicitModel_->get_p_space(p_index),
430 appImplicitModel_->get_f_space()->dim());
431 appImplicitOutArgs.set_DfDp(
432 p_index, MEB::Derivative<Scalar>(my_dfdp_mv_,
433 MEB::DERIV_MV_GRADIENT_FORM));
434 my_dfdp_op_ = my_dfdp_mv_;
435 trans = Thyra::TRANS;
438 TEUCHOS_TEST_FOR_EXCEPTION(
439 true, std::logic_error,
"Invalid df/dp support");
441 appImplicitModel_->evalModel(appImplicitInArgs, appImplicitOutArgs);
444 RCP<DMVPV> f_dfdp = rcp_dynamic_cast<DMVPV>(outArgs.get_f(),
true);
445 const int num_param = f_dfdp->getNonconstMultiVector()->domain()->dim()-1;
446 RCP<Thyra::MultiVectorBase<Scalar> > dfdp =
447 f_dfdp->getNonconstMultiVector()->subView(Range1D(1,num_param));
448 my_dfdp_op_->apply(trans, *dydp, dfdp.ptr(), Scalar(1.0), Scalar(1.0));