Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
Tempus_AdjointSensitivityModelEvaluator_impl.hpp
Go to the documentation of this file.
1// @HEADER
2// ****************************************************************************
3// Tempus: Copyright (2017) Sandia Corporation
4//
5// Distributed under BSD 3-clause license (See accompanying file Copyright.txt)
6// ****************************************************************************
7// @HEADER
8
9#ifndef Tempus_AdjointSensitivityModelEvaluator_impl_hpp
10#define Tempus_AdjointSensitivityModelEvaluator_impl_hpp
11
12#include "Teuchos_TimeMonitor.hpp"
13#include "Tempus_config.hpp"
14
15#include "Thyra_DefaultMultiVectorProductVectorSpace.hpp"
16#include "Thyra_DefaultMultiVectorProductVector.hpp"
19#include "Thyra_DefaultScaledAdjointLinearOp.hpp"
20#include "Thyra_DefaultAdjointLinearOpWithSolve.hpp"
22#include "Thyra_VectorStdOps.hpp"
23#include "Thyra_MultiVectorStdOps.hpp"
24
25namespace Tempus {
26
27template <typename Scalar>
30 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> > & model,
31 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> > & adjoint_residual_model,
32 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> > & adjoint_solve_model,
33 const Scalar& t_init,
34 const Scalar& t_final,
35 const bool is_pseudotransient,
36 const Teuchos::RCP<const Teuchos::ParameterList>& pList) :
37 model_(model),
38 adjoint_residual_model_(adjoint_residual_model),
39 adjoint_solve_model_(adjoint_solve_model),
40 t_init_(t_init),
41 t_final_(t_final),
42 is_pseudotransient_(is_pseudotransient),
43 mass_matrix_is_computed_(false),
44 jacobian_matrix_is_computed_(false),
45 response_gradient_is_computed_(false),
46 t_interp_(Teuchos::ScalarTraits<Scalar>::rmax())
47{
48 typedef Thyra::ModelEvaluatorBase MEB;
49
50 // Set parameters
51 Teuchos::RCP<Teuchos::ParameterList> pl =
52 Teuchos::rcp(new Teuchos::ParameterList);
53 if (pList != Teuchos::null)
54 *pl = *pList;
55 pl->validateParametersAndSetDefaults(*this->getValidParameters());
56 mass_matrix_is_constant_ = pl->get<bool>("Mass Matrix Is Constant");
57 mass_matrix_is_identity_ = pl->get<bool>("Mass Matrix Is Identity");
58 p_index_ = pl->get<int>("Sensitivity Parameter Index", 0);
59 g_index_ = pl->get<int>("Response Function Index", 0);
60 num_adjoint_ = model_->get_g_space(g_index_)->dim();
61
62 // We currently do not support a non-constant mass matrix
63 TEUCHOS_TEST_FOR_EXCEPTION(
64 mass_matrix_is_constant_ == false, std::logic_error,
65 "AdjointSensitivityModelEvaluator currently does not support " <<
66 "non-constant mass matrix df/dx_dot!");
67
69 Thyra::multiVectorProductVectorSpace(model_->get_f_space(), num_adjoint_);
71 Thyra::multiVectorProductVectorSpace(model_->get_x_space(), num_adjoint_);
73 Thyra::multiVectorProductVectorSpace(model_->get_p_space(p_index_),
75
76 // forward and adjoint models must support same InArgs
77 MEB::InArgs<Scalar> me_inArgs = model_->createInArgs();
78 me_inArgs.assertSameSupport(adjoint_residual_model_->createInArgs());
79 me_inArgs.assertSameSupport(adjoint_solve_model_->createInArgs());
80
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 if (me_inArgs.supports(MEB::IN_ARG_alpha))
88 inArgs.setSupports(MEB::IN_ARG_alpha);
89 if (me_inArgs.supports(MEB::IN_ARG_beta))
90 inArgs.setSupports(MEB::IN_ARG_beta);
91
92 // Support additional parameters for x and xdot
93 inArgs.set_Np(me_inArgs.Np());
94 prototypeInArgs_ = inArgs;
95
96 MEB::OutArgs<Scalar> me_outArgs = model_->createOutArgs();
97 MEB::OutArgs<Scalar> adj_mer_outArgs = adjoint_residual_model_->createOutArgs();
98 MEB::OutArgs<Scalar> adj_mes_outArgs = adjoint_solve_model_->createOutArgs();
99 MEB::OutArgsSetup<Scalar> outArgs;
100 outArgs.setModelEvalDescription(this->description());
101 outArgs.set_Np_Ng(me_inArgs.Np(),2);
102 outArgs.setSupports(MEB::OUT_ARG_f);
103 if (adj_mes_outArgs.supports(MEB::OUT_ARG_W_op))
104 outArgs.setSupports(MEB::OUT_ARG_W_op);
105 prototypeOutArgs_ = outArgs;
106
107 // Adjoint residual ME must support W_op to define adjoint ODE/DAE.
108 // Must support alpha, beta if it suports x_dot
109 TEUCHOS_ASSERT(me_inArgs.supports(MEB::IN_ARG_x));
110 TEUCHOS_ASSERT(adj_mer_outArgs.supports(MEB::OUT_ARG_W_op));
111 if (me_inArgs.supports(MEB::IN_ARG_x_dot)) {
112 TEUCHOS_ASSERT(me_inArgs.supports(MEB::IN_ARG_alpha));
113 TEUCHOS_ASSERT(me_inArgs.supports(MEB::IN_ARG_beta));
114 }
115 MEB::DerivativeSupport dgdx_support =
116 me_outArgs.supports(MEB::OUT_ARG_DgDx, g_index_);
117 MEB::DerivativeSupport dgdp_support =
118 me_outArgs.supports(MEB::OUT_ARG_DgDp, g_index_, p_index_);
119 TEUCHOS_ASSERT(dgdx_support.supports(MEB::DERIV_MV_GRADIENT_FORM));
120 TEUCHOS_ASSERT(dgdp_support.supports(MEB::DERIV_MV_GRADIENT_FORM));
121}
122
123template <typename Scalar>
124void
126setFinalTime(const Scalar t_final)
127{
128 t_final_ = t_final;
129}
130
131template <typename Scalar>
132void
135 const Teuchos::RCP<const Tempus::SolutionHistory<Scalar> >& sh)
136{
137 sh_ = sh;
138 if (is_pseudotransient_)
139 forward_state_ = sh_->getCurrentState();
140 else {
141 t_interp_ = Teuchos::ScalarTraits<Scalar>::rmax();
142 forward_state_ = Teuchos::null;
143 }
144
145 // Reset computation flags because we have done a new forward integration
146 mass_matrix_is_computed_ = false;
147 jacobian_matrix_is_computed_ = false;
148 response_gradient_is_computed_ = false;
149}
150
151template <typename Scalar>
152Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
154get_p_space(int p) const
155{
156 TEUCHOS_ASSERT(p < model_->Np());
157 return model_->get_p_space(p);
158}
159
160template <typename Scalar>
161Teuchos::RCP<const Teuchos::Array<std::string> >
163get_p_names(int p) const
164{
165 TEUCHOS_ASSERT(p < model_->Np());
166 return model_->get_p_names(p);
167}
168
169template <typename Scalar>
170Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
172get_x_space() const
173{
174 return adjoint_space_;
175}
176
177template <typename Scalar>
178Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
180get_f_space() const
181{
182 return residual_space_;
183}
184
185template <typename Scalar>
186Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
188get_g_space(int j) const
189{
190 TEUCHOS_ASSERT(j == 0 || j == 1);
191 if (j == 0)
192 return response_space_;
193 return model_->get_g_space(g_index_);
194}
195
196template <typename Scalar>
197Teuchos::RCP<Thyra::LinearOpBase<Scalar> >
199create_W_op() const
200{
201 Teuchos::RCP<Thyra::LinearOpBase<Scalar> > adjoint_op =
202 adjoint_solve_model_->create_W_op();
203 if (adjoint_op == Teuchos::null)
204 return Teuchos::null;
205
206 return Thyra::nonconstMultiVectorLinearOp(adjoint_op, num_adjoint_);
207}
208
209template <typename Scalar>
210Teuchos::RCP<const Thyra::LinearOpWithSolveFactoryBase<Scalar> >
212get_W_factory() const
213{
214 using Teuchos::RCP;
215 using Teuchos::rcp_dynamic_cast;
217
218 RCP<const LOWSFB> alowsfb = adjoint_solve_model_->get_W_factory();
219 if (alowsfb == Teuchos::null)
220 return Teuchos::null; // adjoint_solve_model_ doesn't support W_factory
221
222 return Thyra::multiVectorLinearOpWithSolveFactory(
223 alowsfb, residual_space_, adjoint_space_);
224}
225
226template <typename Scalar>
227Thyra::ModelEvaluatorBase::InArgs<Scalar>
229createInArgs() const
230{
231 return prototypeInArgs_;
232}
233
234template <typename Scalar>
235Thyra::ModelEvaluatorBase::InArgs<Scalar>
237getNominalValues() const
238{
239 typedef Thyra::ModelEvaluatorBase MEB;
240 using Teuchos::RCP;
241 using Teuchos::rcp_dynamic_cast;
242
243 MEB::InArgs<Scalar> me_nominal = model_->getNominalValues();
244 MEB::InArgs<Scalar> nominal = this->createInArgs();
245
246 const Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
247
248 // Set initial x, x_dot
249 RCP< Thyra::VectorBase<Scalar> > x = Thyra::createMember(*adjoint_space_);
250 Thyra::assign(x.ptr(), zero);
251 nominal.set_x(x);
252
253 if (me_nominal.supports(MEB::IN_ARG_x_dot)) {
254 RCP< Thyra::VectorBase<Scalar> > x_dot =
255 Thyra::createMember(*adjoint_space_);
256 Thyra::assign(x_dot.ptr(), zero);
257 nominal.set_x_dot(x_dot);
258 }
259
260 const int np = model_->Np();
261 for (int i=0; i<np; ++i)
262 nominal.set_p(i, me_nominal.get_p(i));
263
264 return nominal;
265}
266
267template <typename Scalar>
268Thyra::ModelEvaluatorBase::OutArgs<Scalar>
270createOutArgsImpl() const
271{
272 return prototypeOutArgs_;
273}
274
275template <typename Scalar>
276void
278evalModelImpl(const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs,
279 const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs) const
280{
281 typedef Thyra::ModelEvaluatorBase MEB;
282 using Teuchos::RCP;
283 using Teuchos::rcp_dynamic_cast;
284
285 TEMPUS_FUNC_TIME_MONITOR_DIFF("Tempus::AdjointSensitivityModelEvaluator::evalModel()", TEMPUS_EVAL);
286
287 // Note: adjoint models compute the transposed W (either explicitly or
288 // implicitly. Thus we need to always call their evalModel() functions
289 // whenever computing the adjoint operators, and subsequent calls to apply()
290 // do not transpose them.
291
292 // Interpolate forward solution at supplied time, reusing previous
293 // interpolation if possible
294 TEUCHOS_ASSERT(sh_ != Teuchos::null);
295 const Scalar t = inArgs.get_t();
296 Scalar forward_t;
297 if (is_pseudotransient_)
298 forward_t = forward_state_->getTime();
299 else {
300 forward_t = t_final_ + t_init_ - t;
301 if (forward_state_ == Teuchos::null || t_interp_ != t) {
302 if (forward_state_ == Teuchos::null)
303 forward_state_ = sh_->interpolateState(forward_t);
304 else
305 sh_->interpolateState(forward_t, forward_state_.get());
306 t_interp_ = t;
307 }
308 }
309
310 // setup input arguments for model
311 MEB::InArgs<Scalar> me_inArgs = model_->getNominalValues();
312 me_inArgs.set_x(forward_state_->getX());
313 if (me_inArgs.supports(MEB::IN_ARG_x_dot)) {
314 if (inArgs.get_x_dot() != Teuchos::null)
315 me_inArgs.set_x_dot(forward_state_->getXDot());
316 else {
317 if (is_pseudotransient_) {
318 // For pseudo-transient, we have to always use the same form of the
319 // residual in order to reuse df/dx, df/dx_dot,..., so we force
320 // the underlying ME to always compute the implicit form with x_dot == 0
321 // if it wasn't provided.
322 if (my_x_dot_ == Teuchos::null) {
323 my_x_dot_ = Thyra::createMember(model_->get_x_space());
324 Thyra::assign(my_x_dot_.ptr(), Scalar(0.0));
325 }
326 me_inArgs.set_x_dot(my_x_dot_);
327 }
328 else {
329 // clear out xdot if it was set in nominalValues to get to ensure we
330 // get the explicit form
331 me_inArgs.set_x_dot(Teuchos::null);
332 }
333 }
334 }
335 if (me_inArgs.supports(MEB::IN_ARG_t))
336 me_inArgs.set_t(forward_t);
337 const int np = me_inArgs.Np();
338 for (int i=0; i<np; ++i)
339 me_inArgs.set_p(i, inArgs.get_p(i));
340
341 // compute adjoint W == model W
342 // It would be nice to not reevaluate W in the psuedo-transient case, but
343 // it isn't clear how to do this in a clean way. Probably just need to
344 // control that with the nonlinear solver.
345 RCP<Thyra::LinearOpBase<Scalar> > op;
346 if (outArgs.supports(MEB::OUT_ARG_W_op))
347 op = outArgs.get_W_op();
348 if (op != Teuchos::null) {
349 TEMPUS_FUNC_TIME_MONITOR_DIFF("Tempus::AdjointSensitivityModelEvaluator::evalModel::W", TEMPUS_EVAL_W);
350 if (me_inArgs.supports(MEB::IN_ARG_alpha))
351 me_inArgs.set_alpha(inArgs.get_alpha());
352 if (me_inArgs.supports(MEB::IN_ARG_beta)) {
353 if (me_inArgs.supports(MEB::IN_ARG_x_dot))
354 me_inArgs.set_beta(inArgs.get_beta()); // Implicit form (see below)
355 else
356 me_inArgs.set_beta(-inArgs.get_beta()); // Explicit form (see below)
357 }
358
359 RCP<Thyra::MultiVectorLinearOp<Scalar> > mv_adjoint_op =
360 rcp_dynamic_cast<Thyra::MultiVectorLinearOp<Scalar> >(op,true);
361 RCP<Thyra::LinearOpBase<Scalar> > adjoint_op =
362 mv_adjoint_op->getNonconstLinearOp();
363 MEB::OutArgs<Scalar> adj_me_outArgs = adjoint_solve_model_->createOutArgs();
364 adj_me_outArgs.set_W_op(adjoint_op);
365 adjoint_solve_model_->evalModel(me_inArgs, adj_me_outArgs);
366 }
367
368 RCP<Thyra::VectorBase<Scalar> > adjoint_f = outArgs.get_f();
369 RCP<Thyra::VectorBase<Scalar> > adjoint_g = outArgs.get_g(0);
370 RCP<Thyra::VectorBase<Scalar> > g = outArgs.get_g(1);
371 RCP<const Thyra::MultiVectorBase<Scalar> > adjoint_x_mv;
372 RCP<const Thyra::VectorBase<Scalar> > adjoint_x;
373 if (adjoint_f != Teuchos::null || adjoint_g != Teuchos::null) {
374 adjoint_x =
375 inArgs.get_x().assert_not_null();
376 adjoint_x_mv =
377 rcp_dynamic_cast<const DMVPV>(adjoint_x,true)->getMultiVector();
378 }
379
380 // Compute adjoint residual F(y):
381 // * For implicit form, F(y) = d/dt( df/dx_dot^T*y ) + df/dx^T*y - dg/dx^T
382 // * For explict form, F(y) = -df/dx^T*y + dg/dx^T
383 // For implicit form, we assume df/dx_dot is constant w.r.t. x, x_dot, and t,
384 // so the residual becomes F(y) = df/dx_dot^T*y_dot + df/dx^T*y - dg/dx^T
385 if (adjoint_f != Teuchos::null) {
386 TEMPUS_FUNC_TIME_MONITOR_DIFF("Tempus::AdjointSensitivityModelEvaluator::evalModel::f", TEMPUS_EVAL_F);
387
388 RCP<Thyra::MultiVectorBase<Scalar> > adjoint_f_mv =
389 rcp_dynamic_cast<DMVPV>(adjoint_f,true)->getNonconstMultiVector();
390
391 MEB::OutArgs<Scalar> me_outArgs = model_->createOutArgs();
392 MEB::OutArgs<Scalar> adj_me_outArgs = adjoint_residual_model_->createOutArgs();
393
394 // dg/dx^T
395 // Don't re-evaluate dg/dx for pseudotransient
396 {
397 TEMPUS_FUNC_TIME_MONITOR_DIFF("Tempus::AdjointSensitivityModelEvaluator::evalModel::dg/dx", TEMPUS_EVAL_DGDX);
398 if (my_dgdx_mv_ == Teuchos::null)
399 my_dgdx_mv_ =
400 Thyra::createMembers(model_->get_x_space(),
401 model_->get_g_space(g_index_)->dim());
402 if (!response_gradient_is_computed_) {
403 me_outArgs.set_DgDx(g_index_,
404 MEB::Derivative<Scalar>(my_dgdx_mv_,
405 MEB::DERIV_MV_GRADIENT_FORM));
406 model_->evalModel(me_inArgs, me_outArgs);
407 me_outArgs.set_DgDx(g_index_, MEB::Derivative<Scalar>());
408 if (is_pseudotransient_)
409 response_gradient_is_computed_ = true;
410 }
411 Thyra::assign(adjoint_f_mv.ptr(), *my_dgdx_mv_);
412 }
413
414 // Explicit form of the residual F(y) = -df/dx^T*y + dg/dx^T
415 // Don't re-evaluate df/dx for pseudotransient
416 {
417 TEMPUS_FUNC_TIME_MONITOR_DIFF("Tempus::AdjointSensitivityModelEvaluator::evalModel::df/dx", TEMPUS_EVAL_DFDX);
418 if (my_dfdx_ == Teuchos::null)
419 my_dfdx_ = adjoint_residual_model_->create_W_op();
420 if (!jacobian_matrix_is_computed_) {
421 adj_me_outArgs.set_W_op(my_dfdx_);
422 if (me_inArgs.supports(MEB::IN_ARG_alpha))
423 me_inArgs.set_alpha(0.0);
424 if (me_inArgs.supports(MEB::IN_ARG_beta))
425 me_inArgs.set_beta(1.0);
426 adjoint_residual_model_->evalModel(me_inArgs, adj_me_outArgs);
427 if (is_pseudotransient_)
428 jacobian_matrix_is_computed_ = true;
429 }
430 my_dfdx_->apply(Thyra::NOTRANS, *adjoint_x_mv, adjoint_f_mv.ptr(),
431 Scalar(-1.0), Scalar(1.0));
432 }
433
434 // Implicit form residual F(y) df/dx_dot^T*y_dot + df/dx^T*y - dg/dx^T
435 // using the second scalar argument to apply() to change the explicit term
436 // above.
437 // Don't re-evaluate df/dx_dot for pseudotransient
438 if (me_inArgs.supports(MEB::IN_ARG_x_dot)) {
439 TEMPUS_FUNC_TIME_MONITOR_DIFF("Tempus::AdjointSensitivityModelEvaluator::evalModel::df/dx_dot", TEMPUS_EVAL_DFDXDOT);
440 RCP<const Thyra::VectorBase<Scalar> > adjoint_x_dot = inArgs.get_x_dot();
441 if (adjoint_x_dot != Teuchos::null) {
442 RCP<const Thyra::MultiVectorBase<Scalar> > adjoint_x_dot_mv =
443 rcp_dynamic_cast<const DMVPV>(adjoint_x_dot,true)->getMultiVector();
444 if (mass_matrix_is_identity_) {
445 // F = -F + y_dot
446 Thyra::V_StVpV(adjoint_f_mv.ptr(), Scalar(-1.0), *adjoint_f_mv,
447 *adjoint_x_dot_mv);
448 }
449 else {
450 if (my_dfdxdot_ == Teuchos::null)
451 my_dfdxdot_ = adjoint_residual_model_->create_W_op();
452 if (!mass_matrix_is_computed_) {
453 adj_me_outArgs.set_W_op(my_dfdxdot_);
454 me_inArgs.set_alpha(1.0);
455 me_inArgs.set_beta(0.0);
456 adjoint_residual_model_->evalModel(me_inArgs, adj_me_outArgs);
457 if (is_pseudotransient_ || mass_matrix_is_constant_)
458 mass_matrix_is_computed_ = true;
459 }
460 my_dfdxdot_->apply(Thyra::NOTRANS, *adjoint_x_dot_mv,
461 adjoint_f_mv.ptr(), Scalar(1.0), Scalar(-1.0));
462 }
463 }
464 }
465 }
466
467 // Compute g = dg/dp^T - df/dp^T*y for computing the model parameter term in
468 // the adjoint sensitivity formula.
469 // We don't add pseudotransient logic here because this part is only
470 // evaluated once in that case anyway.
471 if (adjoint_g != Teuchos::null) {
472 TEMPUS_FUNC_TIME_MONITOR_DIFF("Tempus::AdjointSensitivityModelEvaluator::evalModel::g", TEMPUS_EVAL_G);
473 RCP<Thyra::MultiVectorBase<Scalar> > adjoint_g_mv =
474 rcp_dynamic_cast<DMVPV>(adjoint_g,true)->getNonconstMultiVector();
475
476 MEB::OutArgs<Scalar> me_outArgs = model_->createOutArgs();
477
478 // dg/dp
479 MEB::DerivativeSupport dgdp_support =
480 me_outArgs.supports(MEB::OUT_ARG_DgDp, g_index_, p_index_);
481 if (dgdp_support.supports(MEB::DERIV_MV_GRADIENT_FORM)) {
482 me_outArgs.set_DgDp(g_index_, p_index_,
483 MEB::Derivative<Scalar>(adjoint_g_mv,
484 MEB::DERIV_MV_GRADIENT_FORM));
485 model_->evalModel(me_inArgs, me_outArgs);
486 }
487 else if (dgdp_support.supports(MEB::DERIV_MV_JACOBIAN_FORM)) {
488 const int num_g = model_->get_g_space(g_index_)->dim();
489 const int num_p = model_->get_p_space(p_index_)->dim();
490 RCP<Thyra::MultiVectorBase<Scalar> > dgdp_trans =
491 createMembers(model_->get_g_space(g_index_), num_p);
492 me_outArgs.set_DgDp(g_index_, p_index_,
493 MEB::Derivative<Scalar>(dgdp_trans,
494 MEB::DERIV_MV_JACOBIAN_FORM));
495 model_->evalModel(me_inArgs, me_outArgs);
496 Thyra::DetachedMultiVectorView<Scalar> dgdp_view(*adjoint_g_mv);
497 Thyra::DetachedMultiVectorView<Scalar> dgdp_trans_view(*dgdp_trans);
498 for (int i=0; i<num_p; ++i)
499 for (int j=0; j<num_g; ++j)
500 dgdp_view(i,j) = dgdp_trans_view(j,i);
501 }
502 else
503 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
504 "Invalid dg/dp support");
505 me_outArgs.set_DgDp(g_index_, p_index_, MEB::Derivative<Scalar>());
506
507 // dg/dp - df/dp^T*y
508 MEB::DerivativeSupport dfdp_support =
509 me_outArgs.supports(MEB::OUT_ARG_DfDp, p_index_);
510 Thyra::EOpTransp trans = Thyra::CONJTRANS;
511 if (dfdp_support.supports(MEB::DERIV_LINEAR_OP)) {
512 if (my_dfdp_op_ == Teuchos::null)
513 my_dfdp_op_ = model_->create_DfDp_op(p_index_);
514 me_outArgs.set_DfDp(p_index_, MEB::Derivative<Scalar>(my_dfdp_op_));
515 trans = Thyra::CONJTRANS;
516 }
517 else if (dfdp_support.supports(MEB::DERIV_MV_JACOBIAN_FORM)) {
518 if (my_dfdp_mv_ == Teuchos::null)
519 my_dfdp_mv_ = createMembers(model_->get_f_space(),
520 model_->get_p_space(p_index_)->dim());
521 me_outArgs.set_DfDp(p_index_,
522 MEB::Derivative<Scalar>(my_dfdp_mv_,
523 MEB::DERIV_MV_JACOBIAN_FORM));
524 my_dfdp_op_ = my_dfdp_mv_;
525 trans = Thyra::CONJTRANS;
526 }
527 else if (dfdp_support.supports(MEB::DERIV_MV_GRADIENT_FORM)) {
528 if (my_dfdp_mv_ == Teuchos::null)
529 my_dfdp_mv_ = createMembers(model_->get_p_space(p_index_),
530 model_->get_f_space()->dim());
531 me_outArgs.set_DfDp(p_index_,
532 MEB::Derivative<Scalar>(my_dfdp_mv_,
533 MEB::DERIV_MV_GRADIENT_FORM));
534 my_dfdp_op_ = my_dfdp_mv_;
535 trans = Thyra::CONJ;
536 }
537 else
538 TEUCHOS_TEST_FOR_EXCEPTION(
539 true, std::logic_error, "Invalid df/dp support");
540 model_->evalModel(me_inArgs, me_outArgs);
541 my_dfdp_op_->apply(trans, *adjoint_x_mv, adjoint_g_mv.ptr(),
542 Scalar(-1.0), Scalar(1.0));
543 }
544
545 if (g != Teuchos::null) {
546 MEB::OutArgs<Scalar> me_outArgs = model_->createOutArgs();
547 me_outArgs.set_g(g_index_, g);
548 model_->evalModel(me_inArgs, me_outArgs);
549 }
550}
551
552template<class Scalar>
553Teuchos::RCP<const Teuchos::ParameterList>
556{
557 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
558 pl->set<int>("Sensitivity Parameter Index", 0);
559 pl->set<int>("Response Function Index", 0);
560 pl->set<bool>("Mass Matrix Is Constant", true);
561 pl->set<bool>("Mass Matrix Is Identity", false);
562 return pl;
563}
564
565} // namespace Tempus
566
567#endif
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_p_space(int p) const
Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > adjoint_residual_model_
Thyra::ModelEvaluatorBase::InArgs< Scalar > getNominalValues() const
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_x_space() const
void setForwardSolutionHistory(const Teuchos::RCP< const Tempus::SolutionHistory< Scalar > > &sh)
Set solution history from forward evaluation.
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_f_space() const
Thyra::ModelEvaluatorBase::OutArgs< Scalar > createOutArgsImpl() const
static Teuchos::RCP< const Teuchos::ParameterList > getValidParameters()
Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > model_
AdjointSensitivityModelEvaluator(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &adjoint_residual_model, const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &adjoint_solve_model, const Scalar &t_init, const Scalar &t_final, const bool is_pseudotransient, const Teuchos::RCP< const Teuchos::ParameterList > &pList=Teuchos::null)
Constructor.
Teuchos::RCP< const Thyra::LinearOpWithSolveFactoryBase< Scalar > > get_W_factory() const
Teuchos::RCP< const Teuchos::Array< std::string > > get_p_names(int p) const
Thyra::ModelEvaluatorBase::InArgs< Scalar > createInArgs() const
void evalModelImpl(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs, const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &outArgs) const
Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > adjoint_solve_model_
void setFinalTime(const Scalar t_final)
Set the final time from the forward evaluation.
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_g_space(int j) const
Teuchos::RCP< Thyra::LinearOpBase< Scalar > > create_W_op() const
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...