Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
Tempus_StaggeredForwardSensitivityModelEvaluator_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_StaggeredForwardSensitivityModelEvaluator_impl_hpp
10#define Tempus_StaggeredForwardSensitivityModelEvaluator_impl_hpp
11
12#include "Thyra_DefaultMultiVectorProductVector.hpp"
16#include "Thyra_VectorStdOps.hpp"
17#include "Thyra_MultiVectorStdOps.hpp"
18
19namespace Tempus {
20
21template <typename Scalar>
24 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> > & model,
25 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> > & sens_residual_model,
26 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> > & sens_solve_model,
27 const bool is_pseudotransient,
28 const Teuchos::RCP<const Teuchos::ParameterList>& pList,
29 const Teuchos::RCP<MultiVector>& dxdp_init,
30 const Teuchos::RCP<MultiVector>& dx_dotdp_init,
31 const Teuchos::RCP<MultiVector>& dx_dotdotdp_init) :
32 model_(model),
33 sens_residual_model_(sens_residual_model),
34 sens_solve_model_(sens_solve_model),
35 dxdp_init_(dxdp_init),
36 dx_dotdp_init_(dx_dotdp_init),
37 dx_dotdotdp_init_(dx_dotdotdp_init),
38 p_index_(0),
39 g_index_(-1),
40 x_tangent_index_(1),
41 xdot_tangent_index_(2),
42 xdotdot_tangent_index_(3),
43 use_dfdp_as_tangent_(false),
44 use_dgdp_as_tangent_(false),
45 num_param_(0),
46 num_response_(0),
47 g_offset_(0),
48 is_pseudotransient_(is_pseudotransient),
49 mass_matrix_is_computed_(false),
50 jacobian_matrix_is_computed_(false),
51 acceleration_matrix_is_computed_(false),
52 residual_sensitivity_is_computed_(false),
53 t_interp_(Teuchos::ScalarTraits<Scalar>::rmax())
54{
55 typedef Thyra::ModelEvaluatorBase MEB;
56
57 // Set parameters
58 Teuchos::RCP<Teuchos::ParameterList> pl =
59 Teuchos::rcp(new Teuchos::ParameterList);
60 if (pList != Teuchos::null)
61 *pl = *pList;
62 pl->validateParametersAndSetDefaults(*this->getValidParameters());
63 use_dfdp_as_tangent_ = pl->get<bool>("Use DfDp as Tangent");
64 use_dgdp_as_tangent_ = pl->get<bool>("Use DgDp as Tangent");
65 p_index_ = pl->get<int>("Sensitivity Parameter Index");
66 g_index_ = pl->get<int>("Response Function Index", -1);
67 x_tangent_index_ = pl->get<int>("Sensitivity X Tangent Index");
68 xdot_tangent_index_ = pl->get<int>("Sensitivity X-Dot Tangent Index");
69 xdotdot_tangent_index_ = pl->get<int>("Sensitivity X-Dot-Dot Tangent Index");
70
71 num_param_ = model_->get_p_space(p_index_)->dim();
72 if (g_index_ >= 0) {
73 num_response_ = model_->get_g_space(g_index_)->dim();
74 g_offset_ = 2;
75 }
77 Thyra::multiVectorProductVectorSpace(model_->get_x_space(), num_param_);
79 Thyra::multiVectorProductVectorSpace(model_->get_f_space(), num_param_);
80 if (g_index_ >= 0)
82 Thyra::multiVectorProductVectorSpace(model_->get_g_space(g_index_),
84
85 // forward and sensitivity models must support same InArgs
86 MEB::InArgs<Scalar> me_inArgs = model_->createInArgs();
87 me_inArgs.assertSameSupport(sens_residual_model_->createInArgs());
88 me_inArgs.assertSameSupport(sens_solve_model_->createInArgs());
89
90 MEB::InArgsSetup<Scalar> inArgs;
91 inArgs.setModelEvalDescription(this->description());
92 inArgs.setSupports(MEB::IN_ARG_x);
93 if (me_inArgs.supports(MEB::IN_ARG_x_dot))
94 inArgs.setSupports(MEB::IN_ARG_x_dot);
95 if (me_inArgs.supports(MEB::IN_ARG_t))
96 inArgs.setSupports(MEB::IN_ARG_t);
97 if (me_inArgs.supports(MEB::IN_ARG_alpha))
98 inArgs.setSupports(MEB::IN_ARG_alpha);
99 if (me_inArgs.supports(MEB::IN_ARG_beta))
100 inArgs.setSupports(MEB::IN_ARG_beta);
101 if (me_inArgs.supports(MEB::IN_ARG_W_x_dot_dot_coeff))
102 inArgs.setSupports(MEB::IN_ARG_W_x_dot_dot_coeff);
103
104 // Support additional parameters for x and xdot
105 inArgs.set_Np(me_inArgs.Np());
106 prototypeInArgs_ = inArgs;
107
108 MEB::OutArgs<Scalar> me_outArgs = model_->createOutArgs();
109 MEB::OutArgs<Scalar> sens_mer_outArgs = sens_residual_model_->createOutArgs();
110 MEB::OutArgs<Scalar> sens_mes_outArgs = sens_solve_model_->createOutArgs();
111 MEB::OutArgsSetup<Scalar> outArgs;
112 outArgs.setModelEvalDescription(this->description());
113 outArgs.set_Np_Ng(me_inArgs.Np(),me_outArgs.Ng()+g_offset_);
114 outArgs.setSupports(MEB::OUT_ARG_f);
115 if (sens_mer_outArgs.supports(MEB::OUT_ARG_W_op) ||
116 sens_mes_outArgs.supports(MEB::OUT_ARG_W_op))
117 outArgs.setSupports(MEB::OUT_ARG_W_op);
118
119 // Response 0 is the reduced dg/dp (no sensitivities supported)
120 // Response 1 is the reduced g (no sensitivities supported)
121 for (int j=0; j<me_outArgs.Ng(); ++j) {
122 outArgs.setSupports(MEB::OUT_ARG_DgDx_dot, j+g_offset_,
123 me_outArgs.supports(MEB::OUT_ARG_DgDx_dot, j));
124 outArgs.setSupports(MEB::OUT_ARG_DgDx, j+g_offset_,
125 me_outArgs.supports(MEB::OUT_ARG_DgDx, j));
126 for (int l=0; l<me_outArgs.Np(); ++l) {
127 outArgs.setSupports(MEB::OUT_ARG_DgDp, j+g_offset_, l,
128 me_outArgs.supports(MEB::OUT_ARG_DgDp, j, l));
129 }
130 }
131 prototypeOutArgs_ = outArgs;
132
133 // Sensitivity residual ME must support W_op to define adjoint ODE/DAE.
134 // Must support alpha, beta if it suports x_dot
135 TEUCHOS_ASSERT(me_inArgs.supports(MEB::IN_ARG_x));
137 TEUCHOS_ASSERT(sens_mer_outArgs.supports(MEB::OUT_ARG_W_op));
138 if (me_inArgs.supports(MEB::IN_ARG_x_dot)) {
139 TEUCHOS_ASSERT(me_inArgs.supports(MEB::IN_ARG_alpha));
140 TEUCHOS_ASSERT(me_inArgs.supports(MEB::IN_ARG_beta));
141 }
142 TEUCHOS_ASSERT(me_outArgs.supports(MEB::OUT_ARG_DfDp, p_index_).supports(MEB::DERIV_MV_JACOBIAN_FORM));
143}
144
145template <typename Scalar>
146void
149 const Teuchos::RCP<const Tempus::SolutionHistory<Scalar> >& sh)
150{
151 sh_ = sh;
152 t_interp_ = Teuchos::ScalarTraits<Scalar>::rmax();
153}
154
155template <typename Scalar>
156void
159 const Teuchos::RCP<const Tempus::SolutionState<Scalar> >& s)
160{
161 sh_ = Teuchos::null;
162 forward_state_ = s;
163}
164
165template <typename Scalar>
166Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
168get_p_space(int p) const
169{
170 TEUCHOS_ASSERT(p < model_->Np());
171 return model_->get_p_space(p);
172}
173
174template <typename Scalar>
175Teuchos::RCP<const Teuchos::Array<std::string> >
177get_p_names(int p) const
178{
179 TEUCHOS_ASSERT(p < model_->Np());
180 return model_->get_p_names(p);
181}
182
183template <typename Scalar>
184Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
186get_x_space() const
187{
188 return dxdp_space_;
189}
190
191template <typename Scalar>
192Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
194get_f_space() const
195{
196 return dfdp_space_;
197}
198
199template <typename Scalar>
200Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
202get_g_space(int j) const
203{
204 if (g_index_ >= 0) {
205 if (j == 0)
206 return dgdp_space_;
207 else if (j == 1)
208 return model_->get_g_space(g_index_);
209 }
210 return model_->get_g_space(j-g_offset_);
211}
212
213template <typename Scalar>
214Teuchos::ArrayView<const std::string>
216get_g_names(int j) const
217{
218 if (g_index_ >= 0) {
219 if (j == 0) {
220 Teuchos::Array<std::string> names = model_->get_g_names(g_index_);
221 for (int i=0; i<names.size(); ++i)
222 names[i] = names[i] + "_reduced sensitivity";
223 return names();
224 }
225 else if (j == 1)
226 return model_->get_g_names(g_index_);
227 }
228 return model_->get_g_names(j-g_offset_);
229}
230
231template <typename Scalar>
232Teuchos::RCP<Thyra::LinearOpBase<Scalar> >
234create_W_op() const
235{
236 Teuchos::RCP<Thyra::LinearOpBase<Scalar> > op;
237 if (lo_ != Teuchos::null)
238 op = lo_;
239 else
240 op = sens_solve_model_->create_W_op();
241 return Thyra::nonconstMultiVectorLinearOp(op, num_param_);
242}
243
244template <typename Scalar>
245Teuchos::RCP<Thyra::LinearOpBase<Scalar> >
247create_DgDx_dot_op(int j) const
248{
249 return model_->create_DgDx_dot_op(j-g_offset_);
250}
251
252template <typename Scalar>
253Teuchos::RCP<Thyra::LinearOpBase<Scalar> >
255create_DgDx_op(int j) const
256{
257 return model_->create_DgDx_op(j-g_offset_);
258}
259
260template <typename Scalar>
261Teuchos::RCP<Thyra::LinearOpBase<Scalar> >
263create_DgDp_op(int j, int l) const
264{
265 return model_->create_DgDp_op(j-g_offset_,l);
266}
267
268template <typename Scalar>
269Teuchos::RCP<const Thyra::LinearOpWithSolveFactoryBase<Scalar> >
271get_W_factory() const
272{
273 using Teuchos::RCP;
274 using Teuchos::rcp_dynamic_cast;
275
276 Teuchos::RCP<const Thyra::LinearOpWithSolveFactoryBase<Scalar> > factory;
277 Teuchos::RCP<const Thyra::LinearOpWithSolveFactoryBase<Scalar> > model_factory
278 = sens_solve_model_->get_W_factory();
279 if (model_factory == Teuchos::null)
280 return Teuchos::null; // model_ doesn't support W_factory
281 if (po_ != Teuchos::null) {
282 factory =
283 Thyra::reuseLinearOpWithSolveFactory<Scalar>(model_factory, po_);
284 }
285 else
286 factory = model_factory;
287 return Thyra::multiVectorLinearOpWithSolveFactory(factory,
288 dfdp_space_,
289 dxdp_space_);
290}
291
292template <typename Scalar>
293Thyra::ModelEvaluatorBase::InArgs<Scalar>
295createInArgs() const
296{
297 return prototypeInArgs_;
298}
299
300template <typename Scalar>
301Thyra::ModelEvaluatorBase::InArgs<Scalar>
303getNominalValues() const
304{
305 typedef Thyra::ModelEvaluatorBase MEB;
306 typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;
307 using Teuchos::RCP;
308 using Teuchos::rcp_dynamic_cast;
309
310 MEB::InArgs<Scalar> me_nominal = model_->getNominalValues();
311 MEB::InArgs<Scalar> nominal = this->createInArgs();
312
313 const Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
314
315 // Set initial x. If dxdp_init == null, set initial dx/dp = 0
316 RCP< Thyra::VectorBase<Scalar> > x = Thyra::createMember(*dxdp_space_);
317 RCP<DMVPV> dxdp = rcp_dynamic_cast<DMVPV>(x,true);
318 if (dxdp_init_ == Teuchos::null)
319 Thyra::assign(dxdp->getNonconstMultiVector().ptr(), zero);
320 else
321 Thyra::assign(dxdp->getNonconstMultiVector().ptr(), *dxdp_init_);
322 nominal.set_x(x);
323
324 // Set initial xdot. If dx_dotdp_init == null, set initial dxdot/dp = 0
325 if (me_nominal.supports(MEB::IN_ARG_x_dot)) {
326 RCP< Thyra::VectorBase<Scalar> > xdot = Thyra::createMember(*dxdp_space_);
327 RCP<DMVPV> dxdotdp = rcp_dynamic_cast<DMVPV>(xdot,true);
328 if (dx_dotdp_init_ == Teuchos::null)
329 Thyra::assign(dxdotdp->getNonconstMultiVector().ptr(), zero);
330 else
331 Thyra::assign(dxdotdp->getNonconstMultiVector().ptr(), *dx_dotdp_init_);
332 nominal.set_x_dot(xdot);
333 }
334
335 // Set initial xdotdot. If dx_dotdotdp_init == null, set initial dxdotdot/dp = 0
336 if (me_nominal.supports(MEB::IN_ARG_x_dot_dot)) {
337 RCP< Thyra::VectorBase<Scalar> > xdotdot =
338 Thyra::createMember(*dxdp_space_);
339 RCP<DMVPV> dxdotdotdp = rcp_dynamic_cast<DMVPV>(xdotdot,true);
340 if (dx_dotdotdp_init_ == Teuchos::null)
341 Thyra::assign(dxdotdotdp->getNonconstMultiVector().ptr(), zero);
342 else
343 Thyra::assign(dxdotdotdp->getNonconstMultiVector().ptr(),
344 *dx_dotdotdp_init_);
345 nominal.set_x_dot_dot(xdotdot);
346 }
347
348 const int np = model_->Np();
349 for (int i=0; i<np; ++i)
350 nominal.set_p(i, me_nominal.get_p(i));
351 return nominal;
352}
353
354template <typename Scalar>
355Thyra::ModelEvaluatorBase::OutArgs<Scalar>
357createOutArgsImpl() const
358{
359 return prototypeOutArgs_;
360}
361
362template <typename Scalar>
363void
365evalModelImpl(const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs,
366 const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs) const
367{
368 typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;
369 typedef Thyra::ModelEvaluatorBase MEB;
370 using Teuchos::RCP;
371 using Teuchos::rcp_dynamic_cast;
372
373 // Interpolate forward solution at supplied time, reusing previous
374 // interpolation if possible
375 Scalar forward_t;
376 if (sh_ != Teuchos::null) {
377 forward_t = inArgs.get_t();
378 if (t_interp_ != forward_t) {
379 if (nc_forward_state_ == Teuchos::null)
380 nc_forward_state_ = sh_->interpolateState(forward_t);
381 else
382 sh_->interpolateState(forward_t, nc_forward_state_.get());
383 forward_state_ = nc_forward_state_;
384 t_interp_ = forward_t;
385 }
386 }
387 else {
388 TEUCHOS_ASSERT(forward_state_ != Teuchos::null);
389 forward_t = forward_state_->getTime();
390 }
391
392 const bool use_tangent = use_dfdp_as_tangent_ || use_dgdp_as_tangent_;
393
394 // setup input arguments for model
395 RCP< const Thyra::MultiVectorBase<Scalar> > dxdp, dxdotdp, dxdotdotdp;
396 MEB::InArgs<Scalar> me_inArgs = model_->getNominalValues();
397 dxdp = rcp_dynamic_cast<const DMVPV>(inArgs.get_x(),true)->getMultiVector();
398 me_inArgs.set_x(forward_state_->getX());
399 if (use_dfdp_as_tangent_)
400 me_inArgs.set_p(x_tangent_index_, inArgs.get_x());
401 if (me_inArgs.supports(MEB::IN_ARG_x_dot)) {
402 if (inArgs.get_x_dot() != Teuchos::null) {
403 dxdotdp =
404 rcp_dynamic_cast<const DMVPV>(inArgs.get_x_dot(),true)->getMultiVector();
405 me_inArgs.set_x_dot(forward_state_->getXDot());
406 if (use_dfdp_as_tangent_)
407 me_inArgs.set_p(xdot_tangent_index_, inArgs.get_x_dot());
408 }
409 else {
410 // clear out xdot if it was set in nominalValues to get to ensure we
411 // get the explicit form
412 me_inArgs.set_x_dot(Teuchos::null);
413 }
414 }
415 if (me_inArgs.supports(MEB::IN_ARG_x_dot_dot)) {
416 if (inArgs.get_x_dot_dot() != Teuchos::null) {
417 dxdotdotdp =
418 rcp_dynamic_cast<const DMVPV>(inArgs.get_x_dot_dot(),true)->getMultiVector();
419 me_inArgs.set_x_dot_dot(forward_state_->getXDotDot());
420 if (use_dfdp_as_tangent_)
421 me_inArgs.set_p(xdotdot_tangent_index_, inArgs.get_x_dot_dot());
422 }
423 else // clear out xdotdot if it was set in nominalValues
424 me_inArgs.set_x_dot_dot(Teuchos::null);
425 }
426 if (me_inArgs.supports(MEB::IN_ARG_t))
427 me_inArgs.set_t(forward_t);
428 if (me_inArgs.supports(MEB::IN_ARG_alpha))
429 me_inArgs.set_alpha(inArgs.get_alpha());
430 if (me_inArgs.supports(MEB::IN_ARG_beta))
431 me_inArgs.set_beta(inArgs.get_beta());
432 if (me_inArgs.supports(MEB::IN_ARG_W_x_dot_dot_coeff))
433 me_inArgs.set_W_x_dot_dot_coeff(inArgs.get_W_x_dot_dot_coeff());
434
435 // Set parameters -- be careful to only set ones that were set in our
436 // inArgs to not null out any specified through nominalValues or
437 // dx/dp above
438 const int np = me_inArgs.Np();
439 for (int i=0; i<np; ++i) {
440 if (inArgs.get_p(i) != Teuchos::null)
441 if (!use_tangent ||
442 (use_tangent && i != x_tangent_index_ &&
443 i != xdot_tangent_index_ &&
444 i != xdotdot_tangent_index_ ))
445 me_inArgs.set_p(i, inArgs.get_p(i));
446 }
447
448 // setup output arguments for model
449 RCP< Thyra::MultiVectorBase<Scalar> > dfdp;
450 MEB::OutArgs<Scalar> me_outArgs = model_->createOutArgs();
451 if (outArgs.get_f() != Teuchos::null) {
452 dfdp = rcp_dynamic_cast<DMVPV>(outArgs.get_f(),true)->getNonconstMultiVector();
453 if (!residual_sensitivity_is_computed_) {
454 if (!use_dfdp_as_tangent_ && is_pseudotransient_) {
455 if (my_dfdp_ == Teuchos::null)
456 my_dfdp_ = Thyra::createMembers(model_->get_f_space(),
457 model_->get_p_space(p_index_)->dim());
458 me_outArgs.set_DfDp(p_index_, my_dfdp_);
459 }
460 else
461 me_outArgs.set_DfDp(p_index_, dfdp);
462 }
463 }
464 if (outArgs.supports(MEB::OUT_ARG_W_op) && lo_ == Teuchos::null &&
465 outArgs.get_W_op() != Teuchos::null &&
466 model_.ptr() == sens_solve_model_.ptr()) {
467 RCP<Thyra::LinearOpBase<Scalar> > op = outArgs.get_W_op();
468 RCP<Thyra::MultiVectorLinearOp<Scalar> > mv_op =
469 rcp_dynamic_cast<Thyra::MultiVectorLinearOp<Scalar> >(op,true);
470 me_outArgs.set_W_op(mv_op->getNonconstLinearOp());
471 }
472 for (int j=g_offset_; j<outArgs.Ng(); ++j) {
473 me_outArgs.set_g(j-g_offset_, outArgs.get_g(j));
474 if (!me_outArgs.supports(MEB::OUT_ARG_DgDx_dot,j-g_offset_).none())
475 me_outArgs.set_DgDx_dot(j-g_offset_, outArgs.get_DgDx_dot(j));
476 if (!me_outArgs.supports(MEB::OUT_ARG_DgDx,j-g_offset_).none())
477 me_outArgs.set_DgDx(j-g_offset_, outArgs.get_DgDx(j));
478 for (int l=0; l<outArgs.Np(); ++l)
479 if (!me_outArgs.supports(MEB::OUT_ARG_DgDp,j-g_offset_,l).none())
480 me_outArgs.set_DgDp(j-g_offset_, l, outArgs.get_DgDp(j,l));
481 }
482 if (g_index_ >= 0 && outArgs.get_g(1) != Teuchos::null)
483 me_outArgs.set_g(g_index_, outArgs.get_g(1));
484
485 // build residual and jacobian
486 model_->evalModel(me_inArgs, me_outArgs);
487
488 // Compute W_op separately if we have a separate sensitivity solve ME
489 if (outArgs.supports(MEB::OUT_ARG_W_op) && lo_ == Teuchos::null &&
490 outArgs.get_W_op() != Teuchos::null &&
491 model_.ptr() != sens_solve_model_.ptr()) {
492 MEB::OutArgs<Scalar> sens_me_outArgs = sens_solve_model_->createOutArgs();
493 RCP<Thyra::LinearOpBase<Scalar> > op = outArgs.get_W_op();
494 RCP<Thyra::MultiVectorLinearOp<Scalar> > mv_op =
495 rcp_dynamic_cast<Thyra::MultiVectorLinearOp<Scalar> >(op,true);
496 sens_me_outArgs.set_W_op(mv_op->getNonconstLinearOp());
497 sens_solve_model_->evalModel(me_inArgs, sens_me_outArgs);
498 }
499
500 // Compute (df/dx) * (dx/dp) + (df/dxdot) * (dxdot/dp) + (df/dxdotdot) * (dxdotdot/dp) + (df/dp)
501 // if the underlying ME doesn't already do this. This requires computing
502 // df/dx, df/dxdot, df/dxdotdot as separate operators.
503 // For pseudo-transient, we would like to reuse these operators, but this is
504 // complicated when steppers use both implicit and explicit forms.
505 if (!use_dfdp_as_tangent_) {
506 if (dfdp != Teuchos::null && is_pseudotransient_) {
507 residual_sensitivity_is_computed_ = true;
508 Thyra::assign(dfdp.ptr(), *my_dfdp_);
509 }
510 if (dxdp != Teuchos::null && dfdp != Teuchos::null) {
511 if (my_dfdx_ == Teuchos::null)
512 my_dfdx_ = sens_residual_model_->create_W_op();
513 if (!jacobian_matrix_is_computed_) {
514 MEB::OutArgs<Scalar> meo = sens_residual_model_->createOutArgs();
515 meo.set_W_op(my_dfdx_);
516 if (me_inArgs.supports(MEB::IN_ARG_alpha))
517 me_inArgs.set_alpha(0.0);
518 if (me_inArgs.supports(MEB::IN_ARG_beta))
519 me_inArgs.set_beta(1.0);
520 if (me_inArgs.supports(MEB::IN_ARG_W_x_dot_dot_coeff))
521 me_inArgs.set_W_x_dot_dot_coeff(0.0);
522 sens_residual_model_->evalModel(me_inArgs, meo);
523 if (is_pseudotransient_)
524 jacobian_matrix_is_computed_ = true;
525 }
526 my_dfdx_->apply(Thyra::NOTRANS, *dxdp, dfdp.ptr(), Scalar(1.0), Scalar(1.0));
527 }
528 if (dxdotdp != Teuchos::null && dfdp != Teuchos::null) {
529 if (my_dfdxdot_ == Teuchos::null)
530 my_dfdxdot_ = sens_residual_model_->create_W_op();
531 if (!mass_matrix_is_computed_) {
532 MEB::OutArgs<Scalar> meo = sens_residual_model_->createOutArgs();
533 meo.set_W_op(my_dfdxdot_);
534 if (me_inArgs.supports(MEB::IN_ARG_alpha))
535 me_inArgs.set_alpha(1.0);
536 if (me_inArgs.supports(MEB::IN_ARG_beta))
537 me_inArgs.set_beta(0.0);
538 if (me_inArgs.supports(MEB::IN_ARG_W_x_dot_dot_coeff))
539 me_inArgs.set_W_x_dot_dot_coeff(0.0);
540 sens_residual_model_->evalModel(me_inArgs, meo);
541 if (is_pseudotransient_)
542 mass_matrix_is_computed_ = true;
543 }
544 my_dfdxdot_->apply(Thyra::NOTRANS, *dxdotdp, dfdp.ptr(), Scalar(1.0), Scalar(1.0));
545 }
546 if (dxdotdotdp != Teuchos::null && dfdp != Teuchos::null) {
547 if (my_dfdxdotdot_ == Teuchos::null)
548 my_dfdxdotdot_ = sens_residual_model_->create_W_op();
549 if (!acceleration_matrix_is_computed_) {
550 MEB::OutArgs<Scalar> meo = sens_residual_model_->createOutArgs();
551 meo.set_W_op(my_dfdxdotdot_);
552 if (me_inArgs.supports(MEB::IN_ARG_alpha))
553 me_inArgs.set_alpha(0.0);
554 if (me_inArgs.supports(MEB::IN_ARG_beta))
555 me_inArgs.set_beta(0.0);
556 if (me_inArgs.supports(MEB::IN_ARG_W_x_dot_dot_coeff))
557 me_inArgs.set_W_x_dot_dot_coeff(1.0);
558 sens_residual_model_->evalModel(me_inArgs, meo);
559 if (is_pseudotransient_)
560 acceleration_matrix_is_computed_ = true;
561 }
562 my_dfdxdotdot_->apply(Thyra::NOTRANS, *dxdotdotdp, dfdp.ptr(), Scalar(1.0), Scalar(1.0));
563 }
564 }
565
566 if (g_index_ >= 0 && outArgs.get_g(0) != Teuchos::null) {
567 MEB::OutArgs<Scalar> meo = model_->createOutArgs();
568 RCP<Thyra::MultiVectorBase<Scalar> > dgdp =
569 rcp_dynamic_cast<DMVPV>(outArgs.get_g(0),true)->getNonconstMultiVector();
570 RCP<Thyra::MultiVectorBase<Scalar> > dgdp_trans;
571 MEB::DerivativeSupport dgdp_support =
572 meo.supports(MEB::OUT_ARG_DgDp, g_index_, p_index_);
573 if (dgdp_support.supports(MEB::DERIV_MV_JACOBIAN_FORM))
574 meo.set_DgDp(g_index_, p_index_,
575 MEB::Derivative<Scalar>(dgdp, MEB::DERIV_MV_JACOBIAN_FORM));
576 else if (dgdp_support.supports(MEB::DERIV_MV_GRADIENT_FORM)) {
577 dgdp_trans = createMembers(model_->get_p_space(p_index_),
578 num_response_);
579 meo.set_DgDp(g_index_, p_index_,
580 MEB::Derivative<Scalar>(dgdp_trans, MEB::DERIV_MV_GRADIENT_FORM));
581 }
582 else
583 TEUCHOS_TEST_FOR_EXCEPTION(
584 true, std::logic_error,
585 "Operator form of dg/dp not supported for reduced sensitivity");
586
587 if (!use_dgdp_as_tangent_ || dxdp == Teuchos::null) {
588 MEB::DerivativeSupport dgdx_support =
589 meo.supports(MEB::OUT_ARG_DgDx, g_index_);
590 if (dgdx_support.supports(MEB::DERIV_LINEAR_OP)) {
591 if (my_dgdx_ == Teuchos::null)
592 my_dgdx_ = model_->create_DgDx_op(g_index_);
593 meo.set_DgDx(g_index_, my_dgdx_);
594 }
595 else if (dgdx_support.supports(MEB::DERIV_MV_GRADIENT_FORM)) {
596 if (my_dgdx_mv_ == Teuchos::null)
597 my_dgdx_mv_ = createMembers(model_->get_g_space(g_index_), num_param_);
598 meo.set_DgDx(g_index_,
599 MEB::Derivative<Scalar>(my_dgdx_mv_,
600 MEB::DERIV_MV_GRADIENT_FORM));
601 }
602 else
603 TEUCHOS_TEST_FOR_EXCEPTION(
604 true, std::logic_error,
605 "Jacobian form of dg/dx not supported for reduced sensitivity");
606 }
607
608 // Clear dx/dp, dxdot/dp, dxdotdot/dp from inArgs if set and we are not
609 // using dg/dp as the tangent. Note, even though dxdp, dxdotdp, and
610 // dxdotdotdp are not used here, they are non-null if x/xdot/xdotdot are
611 // supported and supplied.
612 if (!use_dgdp_as_tangent_ && use_dfdp_as_tangent_ &&
613 dxdp != Teuchos::null)
614 me_inArgs.set_p(x_tangent_index_, Teuchos::null);
615 if (!use_dgdp_as_tangent_ && use_dfdp_as_tangent_ &&
616 dxdotdp != Teuchos::null)
617 me_inArgs.set_p(xdot_tangent_index_, Teuchos::null);
618 if (!use_dgdp_as_tangent_ && use_dfdp_as_tangent_ &&
619 dxdotdotdp != Teuchos::null)
620 me_inArgs.set_p(xdotdot_tangent_index_, Teuchos::null);
621
622 // Set dx/dp, dxdot/dp, dxdodot/dp if necessary
623 if (use_dgdp_as_tangent_ && !use_dfdp_as_tangent_ &&
624 dxdp != Teuchos::null)
625 me_inArgs.set_p(x_tangent_index_, inArgs.get_x());
626 if (use_dgdp_as_tangent_ && !use_dfdp_as_tangent_ &&
627 dxdotdp != Teuchos::null)
628 me_inArgs.set_p(xdot_tangent_index_, inArgs.get_x_dot());
629 if (use_dgdp_as_tangent_ && !use_dfdp_as_tangent_ &&
630 dxdotdotdp != Teuchos::null)
631 me_inArgs.set_p(xdotdot_tangent_index_, inArgs.get_x_dot_dot());
632
633 model_->evalModel(me_inArgs, meo);
634
635 // transpose reduced dg/dp if necessary
636 if (dgdp_trans != Teuchos::null) {
637 Thyra::DetachedMultiVectorView<Scalar> dgdp_view(*dgdp);
638 Thyra::DetachedMultiVectorView<Scalar> dgdp_trans_view(*dgdp_trans);
639 for (int i=0; i<num_param_; ++i)
640 for (int j=0; j<num_response_; ++j)
641 dgdp_view(j,i) = dgdp_trans_view(i,j);
642 }
643
644 // Compute (dg/dx) * (dx/dp) + (dg/dp) if the underlying ME doesn't already
645 // do this.
646 if (!use_dgdp_as_tangent_ && dxdp != Teuchos::null) {
647 MEB::DerivativeSupport dgdx_support =
648 me_outArgs.supports(MEB::OUT_ARG_DgDx, g_index_);
649 if (dgdx_support.supports(MEB::DERIV_LINEAR_OP)) {
650 my_dgdx_->apply(Thyra::NOTRANS, *dxdp, dgdp.ptr(), Scalar(1.0), Scalar(1.0));
651 }
652 else if (dgdx_support.supports(MEB::DERIV_MV_GRADIENT_FORM)) {
653 my_dgdx_mv_->apply(Thyra::TRANS, *dxdp, dgdp.ptr(), Scalar(1.0), Scalar(1.0));
654 }
655 else
656 TEUCHOS_TEST_FOR_EXCEPTION(
657 true, std::logic_error,
658 "Jacobian form of dg/dx not supported for reduced sensitivity");
659 }
660 }
661}
662
663template<class Scalar>
664Teuchos::RCP<const Teuchos::ParameterList>
667{
668 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
669 pl->set<bool>("Use DfDp as Tangent", false);
670 pl->set<bool>("Use DgDp as Tangent", false);
671 pl->set<int>("Sensitivity Parameter Index", 0);
672 pl->set<int>("Response Function Index", -1);
673 pl->set<int>("Sensitivity X Tangent Index", 1);
674 pl->set<int>("Sensitivity X-Dot Tangent Index", 2);
675 pl->set<int>("Sensitivity X-Dot-Dot Tangent Index", 3);
676 return pl;
677}
678
679} // namespace Tempus
680
681#endif
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
Solution state for integrators and steppers. SolutionState contains the metadata for solutions and th...
Teuchos::RCP< const Thyra::LinearOpWithSolveFactoryBase< Scalar > > get_W_factory() const
Teuchos::RCP< Thyra::LinearOpBase< Scalar > > create_DgDx_op(int j) const
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_x_space() const
Teuchos::RCP< Thyra::LinearOpBase< Scalar > > create_DgDx_dot_op(int j) const
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_p_space(int p) const
StaggeredForwardSensitivityModelEvaluator(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &sens_residual_model, const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &sens_solve_model, const bool is_pseudotransient, const Teuchos::RCP< const Teuchos::ParameterList > &pList=Teuchos::null, const Teuchos::RCP< MultiVector > &dxdp_init=Teuchos::null, const Teuchos::RCP< MultiVector > &dx_dotdp_init=Teuchos::null, const Teuchos::RCP< MultiVector > &dx_dotdot_dp_init=Teuchos::null)
Constructor.
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_f_space() const
Teuchos::RCP< Thyra::LinearOpBase< Scalar > > create_DgDp_op(int j, int l) const
void setForwardSolutionHistory(const Teuchos::RCP< const Tempus::SolutionHistory< Scalar > > &sh)
Set solution history from forward state evaluation (for interpolation)
void evalModelImpl(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs, const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &outArgs) const
Teuchos::RCP< const Teuchos::Array< std::string > > get_p_names(int p) const
virtual void setForwardSolutionState(const Teuchos::RCP< const Tempus::SolutionState< Scalar > > &s)
Set solution state from forward state evaluation (for frozen state)
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_g_space(int j) const