29 const Scalar& t_final,
30 const Teuchos::RCP<const Teuchos::ParameterList>& pList) :
32 adjoint_model_(adjoint_model),
35 mass_matrix_is_computed_(false),
36 t_interp_(
Teuchos::ScalarTraits<Scalar>::rmax())
40 using Thyra::VectorSpaceBase;
41 typedef Thyra::ModelEvaluatorBase MEB;
44 Teuchos::RCP<Teuchos::ParameterList> pl =
45 Teuchos::rcp(
new Teuchos::ParameterList);
46 if (pList != Teuchos::null)
51 p_index_ = pl->get<
int>(
"Sensitivity Parameter Index", 0);
52 g_index_ = pl->get<
int>(
"Response Function Index", 0);
56 TEUCHOS_TEST_FOR_EXCEPTION(
58 "AdjointAuxSensitivityModelEvaluator currently does not support " <<
59 "non-constant mass matrix df/dx_dot!");
66 Thyra::multiVectorProductVectorSpace(
model_->get_p_space(
p_index_),
68 Array< RCP<const VectorSpaceBase<Scalar> > > x_spaces(2);
69 Array< RCP<const VectorSpaceBase<Scalar> > > f_spaces(2);
78 MEB::InArgs<Scalar> me_inArgs =
model_->createInArgs();
81 MEB::InArgsSetup<Scalar> inArgs;
82 inArgs.setModelEvalDescription(this->description());
83 inArgs.setSupports(MEB::IN_ARG_x);
84 inArgs.setSupports(MEB::IN_ARG_t);
85 if (me_inArgs.supports(MEB::IN_ARG_x_dot))
86 inArgs.setSupports(MEB::IN_ARG_x_dot);
87 inArgs.setSupports(MEB::IN_ARG_alpha);
88 inArgs.setSupports(MEB::IN_ARG_beta);
91 inArgs.set_Np(me_inArgs.Np());
94 MEB::OutArgs<Scalar> me_outArgs =
model_->createOutArgs();
95 MEB::OutArgs<Scalar> adj_me_outArgs =
adjoint_model_->createOutArgs();
96 MEB::OutArgsSetup<Scalar> outArgs;
97 outArgs.setModelEvalDescription(this->description());
98 outArgs.set_Np_Ng(me_inArgs.Np(),0);
99 outArgs.setSupports(MEB::OUT_ARG_f);
100 outArgs.setSupports(MEB::OUT_ARG_W_op);
105 TEUCHOS_ASSERT(me_inArgs.supports(MEB::IN_ARG_x));
106 TEUCHOS_ASSERT(adj_me_outArgs.supports(MEB::OUT_ARG_W_op));
107 if (me_inArgs.supports(MEB::IN_ARG_x_dot)) {
108 TEUCHOS_ASSERT(me_inArgs.supports(MEB::IN_ARG_alpha));
109 TEUCHOS_ASSERT(me_inArgs.supports(MEB::IN_ARG_beta));
263evalModelImpl(
const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs,
264 const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs)
const
266 typedef Thyra::ModelEvaluatorBase MEB;
268 using Teuchos::rcp_dynamic_cast;
277 TEUCHOS_ASSERT(sh_ != Teuchos::null);
278 const Scalar t = inArgs.get_t();
279 const Scalar forward_t = t_final_ + t_init_ - t;
280 if (forward_state_ == Teuchos::null || t_interp_ != t) {
281 if (forward_state_ == Teuchos::null)
282 forward_state_ = sh_->interpolateState(forward_t);
284 sh_->interpolateState(forward_t, forward_state_.get());
289 MEB::InArgs<Scalar> me_inArgs = model_->getNominalValues();
290 me_inArgs.set_x(forward_state_->getX());
291 if (me_inArgs.supports(MEB::IN_ARG_x_dot) &&
292 inArgs.get_x_dot() != Teuchos::null)
293 me_inArgs.set_x_dot(forward_state_->getXDot());
294 if (me_inArgs.supports(MEB::IN_ARG_t))
295 me_inArgs.set_t(forward_t);
296 const int np = me_inArgs.Np();
297 for (
int i=0; i<np; ++i)
298 me_inArgs.set_p(i, inArgs.get_p(i));
301 RCP<Thyra::LinearOpBase<Scalar> > op = outArgs.get_W_op();
302 if (op != Teuchos::null) {
303 if (me_inArgs.supports(MEB::IN_ARG_alpha))
304 me_inArgs.set_alpha(inArgs.get_alpha());
305 if (me_inArgs.supports(MEB::IN_ARG_beta))
306 me_inArgs.set_beta(inArgs.get_beta());
309 RCP<Thyra::DefaultBlockedLinearOp<Scalar> > block_op =
310 rcp_dynamic_cast<Thyra::DefaultBlockedLinearOp<Scalar> >(op,
true);
311 RCP<Thyra::MultiVectorLinearOp<Scalar> > mv_adjoint_op =
312 rcp_dynamic_cast<Thyra::MultiVectorLinearOp<Scalar> >(
313 block_op->getNonconstBlock(0,0),
true);
314 RCP<Thyra::LinearOpBase<Scalar> > adjoint_op =
315 mv_adjoint_op->getNonconstLinearOp();
316 MEB::OutArgs<Scalar> adj_me_outArgs = adjoint_model_->createOutArgs();
317 adj_me_outArgs.set_W_op(adjoint_op);
318 adjoint_model_->evalModel(me_inArgs, adj_me_outArgs);
321 RCP<Thyra::ScaledIdentityLinearOpWithSolve<Scalar> > si_op =
322 rcp_dynamic_cast<Thyra::ScaledIdentityLinearOpWithSolve<Scalar> >(
323 block_op->getNonconstBlock(1,1),
true);
324 si_op->setScale(inArgs.get_alpha());
332 RCP<Thyra::VectorBase<Scalar> > f = outArgs.get_f();
333 if (f != Teuchos::null) {
334 RCP<const Thyra::VectorBase<Scalar> > x = inArgs.get_x().assert_not_null();
335 RCP<const DPV> prod_x = rcp_dynamic_cast<const DPV>(x,
true);
336 RCP<const Thyra::VectorBase<Scalar> > adjoint_x = prod_x->getVectorBlock(0);
337 RCP<const Thyra::MultiVectorBase<Scalar> >adjoint_x_mv =
338 rcp_dynamic_cast<const DMVPV>(adjoint_x,
true)->getMultiVector();
340 RCP<DPV> prod_f = rcp_dynamic_cast<DPV>(f,
true);
341 RCP<Thyra::VectorBase<Scalar> > adjoint_f =
342 prod_f->getNonconstVectorBlock(0);
343 RCP<Thyra::MultiVectorBase<Scalar> > adjoint_f_mv =
344 rcp_dynamic_cast<DMVPV>(adjoint_f,
true)->getNonconstMultiVector();
346 MEB::OutArgs<Scalar> adj_me_outArgs = adjoint_model_->createOutArgs();
348 if (my_dfdx_ == Teuchos::null)
349 my_dfdx_ = adjoint_model_->create_W_op();
350 adj_me_outArgs.set_W_op(my_dfdx_);
351 if (me_inArgs.supports(MEB::IN_ARG_alpha))
352 me_inArgs.set_alpha(0.0);
353 if (me_inArgs.supports(MEB::IN_ARG_beta))
354 me_inArgs.set_beta(1.0);
355 adjoint_model_->evalModel(me_inArgs, adj_me_outArgs);
358 my_dfdx_->apply(Thyra::NOTRANS, *adjoint_x_mv, adjoint_f_mv.ptr(),
359 Scalar(-1.0), Scalar(0.0));
363 RCP<const DPV> prod_x_dot;
364 if (me_inArgs.supports(MEB::IN_ARG_x_dot)) {
365 RCP<const Thyra::VectorBase<Scalar> > x_dot = inArgs.get_x_dot();
366 if (x_dot != Teuchos::null) {
367 prod_x_dot = rcp_dynamic_cast<const DPV>(x_dot,
true);
368 RCP<const Thyra::VectorBase<Scalar> > adjoint_x_dot =
369 prod_x_dot->getVectorBlock(0);
370 RCP<const Thyra::MultiVectorBase<Scalar> > adjoint_x_dot_mv =
371 rcp_dynamic_cast<const DMVPV>(adjoint_x_dot,
true)->getMultiVector();
372 if (mass_matrix_is_identity_) {
374 Thyra::V_StVpV(adjoint_f_mv.ptr(), Scalar(-1.0), *adjoint_f_mv,
378 if (my_dfdxdot_ == Teuchos::null)
379 my_dfdxdot_ = adjoint_model_->create_W_op();
380 if (!mass_matrix_is_constant_ || !mass_matrix_is_computed_) {
381 adj_me_outArgs.set_W_op(my_dfdxdot_);
382 me_inArgs.set_alpha(1.0);
383 me_inArgs.set_beta(0.0);
384 adjoint_model_->evalModel(me_inArgs, adj_me_outArgs);
386 mass_matrix_is_computed_ =
true;
388 my_dfdxdot_->apply(Thyra::NOTRANS, *adjoint_x_dot_mv,
389 adjoint_f_mv.ptr(), Scalar(1.0), Scalar(-1.0));
396 RCP<Thyra::VectorBase<Scalar> > adjoint_g =
397 prod_f->getNonconstVectorBlock(1);
398 RCP<Thyra::MultiVectorBase<Scalar> > adjoint_g_mv =
399 rcp_dynamic_cast<DMVPV>(adjoint_g,
true)->getNonconstMultiVector();
401 MEB::OutArgs<Scalar> me_outArgs2 = model_->createOutArgs();
402 MEB::DerivativeSupport dfdp_support =
403 me_outArgs2.supports(MEB::OUT_ARG_DfDp, p_index_);
404 Thyra::EOpTransp trans = Thyra::CONJTRANS;
405 if (dfdp_support.supports(MEB::DERIV_LINEAR_OP)) {
406 if (my_dfdp_op_ == Teuchos::null)
407 my_dfdp_op_ = model_->create_DfDp_op(p_index_);
408 me_outArgs2.set_DfDp(p_index_, MEB::Derivative<Scalar>(my_dfdp_op_));
409 trans = Thyra::CONJTRANS;
411 else if (dfdp_support.supports(MEB::DERIV_MV_JACOBIAN_FORM)) {
412 if (my_dfdp_mv_ == Teuchos::null)
413 my_dfdp_mv_ = createMembers(model_->get_f_space(),
414 model_->get_p_space(p_index_)->dim());
415 me_outArgs2.set_DfDp(p_index_,
416 MEB::Derivative<Scalar>(my_dfdp_mv_,
417 MEB::DERIV_MV_JACOBIAN_FORM));
418 my_dfdp_op_ = my_dfdp_mv_;
419 trans = Thyra::CONJTRANS;
421 else if (dfdp_support.supports(MEB::DERIV_MV_GRADIENT_FORM)) {
422 if (my_dfdp_mv_ == Teuchos::null)
423 my_dfdp_mv_ = createMembers(model_->get_p_space(p_index_),
424 model_->get_f_space()->dim());
425 me_outArgs2.set_DfDp(p_index_,
426 MEB::Derivative<Scalar>(my_dfdp_mv_,
427 MEB::DERIV_MV_GRADIENT_FORM));
428 my_dfdp_op_ = my_dfdp_mv_;
432 TEUCHOS_TEST_FOR_EXCEPTION(
433 true, std::logic_error,
"Invalid df/dp support");
434 model_->evalModel(me_inArgs, me_outArgs2);
435 my_dfdp_op_->apply(trans, *adjoint_x_mv, adjoint_g_mv.ptr(),
436 Scalar(1.0), Scalar(0.0));
438 if (prod_x_dot != Teuchos::null) {
439 RCP<const Thyra::VectorBase<Scalar> > z_dot =
440 prod_x_dot->getVectorBlock(1);
441 RCP<const Thyra::MultiVectorBase<Scalar> > z_dot_mv =
442 rcp_dynamic_cast<const DMVPV>(z_dot,
true)->getMultiVector();
443 Thyra::V_VmV(adjoint_g_mv.ptr(), *z_dot_mv, *adjoint_g_mv);