Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
Tempus_CombinedForwardSensitivityModelEvaluator_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_CombinedForwardSensitivityModelEvaluator_impl_hpp
10#define Tempus_CombinedForwardSensitivityModelEvaluator_impl_hpp
11
12#include "Thyra_DefaultMultiVectorProductVector.hpp"
13#include "Thyra_MultiVectorLinearOp.hpp"
14#include "Thyra_MultiVectorLinearOpWithSolveFactory.hpp"
15#include "Thyra_ReuseLinearOpWithSolveFactory.hpp"
16#include "Thyra_VectorStdOps.hpp"
17#include "Thyra_MultiVectorStdOps.hpp"
18
19namespace Tempus {
20
21template <typename Scalar>
22CombinedForwardSensitivityModelEvaluator<Scalar>::
23CombinedForwardSensitivityModelEvaluator(
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 Teuchos::RCP<const Teuchos::ParameterList>& pList,
28 const Teuchos::RCP<MultiVector>& dxdp_init,
29 const Teuchos::RCP<MultiVector>& dx_dotdp_init,
30 const Teuchos::RCP<MultiVector>& dx_dotdotdp_init) :
31 model_(model),
32 sens_residual_model_(sens_residual_model),
33 sens_solve_model_(sens_solve_model),
34 dxdp_init_(dxdp_init),
35 dx_dotdp_init_(dx_dotdp_init),
36 dx_dotdotdp_init_(dx_dotdotdp_init),
37 p_index_(0),
38 g_index_(-1),
39 x_tangent_index_(1),
40 xdot_tangent_index_(2),
41 xdotdot_tangent_index_(3),
42 use_dfdp_as_tangent_(false),
43 use_dgdp_as_tangent_(false),
44 num_param_(0),
45 num_response_(0),
46 g_offset_(0)
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 use_dfdp_as_tangent_ = pl->get<bool>("Use DfDp as Tangent");
57 use_dgdp_as_tangent_ = pl->get<bool>("Use DgDp as Tangent");
58 p_index_ = pl->get<int>("Sensitivity Parameter Index");
59 g_index_ = pl->get<int>("Response Function Index");
60 x_tangent_index_ = pl->get<int>("Sensitivity X Tangent Index");
61 xdot_tangent_index_ = pl->get<int>("Sensitivity X-Dot Tangent Index");
62 xdotdot_tangent_index_ = pl->get<int>("Sensitivity X-Dot-Dot Tangent Index");
63
64 num_param_ = model_->get_p_space(p_index_)->dim();
65 if (g_index_ >= 0) {
66 num_response_ = model_->get_g_space(g_index_)->dim();
67 g_offset_ = 2;
68 }
69 dxdp_space_ =
70 Thyra::multiVectorProductVectorSpace(model_->get_x_space(), num_param_);
71 x_dxdp_space_ =
72 Thyra::multiVectorProductVectorSpace(model_->get_x_space(), num_param_+1);
73 dfdp_space_ =
74 Thyra::multiVectorProductVectorSpace(model_->get_f_space(), num_param_);
75 f_dfdp_space_ =
76 Thyra::multiVectorProductVectorSpace(model_->get_f_space(), num_param_+1);
77 if (g_index_ >= 0)
78 dgdp_space_ =
79 Thyra::multiVectorProductVectorSpace(model_->get_g_space(g_index_),
80 num_param_);
81
82 // forward and sensitivity models must support same InArgs
83 MEB::InArgs<Scalar> me_inArgs = model_->createInArgs();
84 me_inArgs.assertSameSupport(sens_residual_model_->createInArgs());
85 me_inArgs.assertSameSupport(sens_solve_model_->createInArgs());
86
87 MEB::InArgsSetup<Scalar> inArgs;
88 inArgs.setModelEvalDescription(this->description());
89 inArgs.setSupports(MEB::IN_ARG_x);
90 if (me_inArgs.supports(MEB::IN_ARG_x_dot))
91 inArgs.setSupports(MEB::IN_ARG_x_dot);
92 if (me_inArgs.supports(MEB::IN_ARG_t))
93 inArgs.setSupports(MEB::IN_ARG_t);
94 if (me_inArgs.supports(MEB::IN_ARG_alpha))
95 inArgs.setSupports(MEB::IN_ARG_alpha);
96 if (me_inArgs.supports(MEB::IN_ARG_beta))
97 inArgs.setSupports(MEB::IN_ARG_beta);
98 if (me_inArgs.supports(MEB::IN_ARG_W_x_dot_dot_coeff))
99 inArgs.setSupports(MEB::IN_ARG_W_x_dot_dot_coeff);
100
101 // Support additional parameters for x and xdot
102 inArgs.set_Np(me_inArgs.Np());
103 prototypeInArgs_ = inArgs;
104
105 MEB::OutArgs<Scalar> me_outArgs = model_->createOutArgs();
106 MEB::OutArgs<Scalar> sens_mer_outArgs = sens_residual_model_->createOutArgs();
107 MEB::OutArgs<Scalar> sens_mes_outArgs = sens_solve_model_->createOutArgs();
108 MEB::OutArgsSetup<Scalar> outArgs;
109 outArgs.setModelEvalDescription(this->description());
110 outArgs.set_Np_Ng(me_outArgs.Np(),me_outArgs.Ng()+g_offset_);
111 outArgs.setSupports(MEB::OUT_ARG_f);
112 if (sens_mer_outArgs.supports(MEB::OUT_ARG_W_op) ||
113 sens_mes_outArgs.supports(MEB::OUT_ARG_W_op))
114 outArgs.setSupports(MEB::OUT_ARG_W_op);
115
116 // Response 0 is the reduced dg/dp (no sensitivities supported)
117 // Response 1 is the reduced g (no sensitivities supported)
118 for (int j=0; j<me_outArgs.Ng(); ++j) {
119 outArgs.setSupports(MEB::OUT_ARG_DgDx_dot, j+g_offset_,
120 me_outArgs.supports(MEB::OUT_ARG_DgDx_dot, j));
121 outArgs.setSupports(MEB::OUT_ARG_DgDx, j+g_offset_,
122 me_outArgs.supports(MEB::OUT_ARG_DgDx, j));
123 for (int l=0; l<me_outArgs.Np(); ++l) {
124 outArgs.setSupports(MEB::OUT_ARG_DgDp, j+g_offset_, l,
125 me_outArgs.supports(MEB::OUT_ARG_DgDp, j, l));
126 }
127 }
128 prototypeOutArgs_ = outArgs;
129
130 // Sensitivity residual ME must support W_op to define adjoint ODE/DAE.
131 // Must support alpha, beta if it suports x_dot
132 TEUCHOS_ASSERT(me_inArgs.supports(MEB::IN_ARG_x));
133 if (!use_dfdp_as_tangent_)
134 TEUCHOS_ASSERT(sens_mer_outArgs.supports(MEB::OUT_ARG_W_op));
135 if (me_inArgs.supports(MEB::IN_ARG_x_dot)) {
136 TEUCHOS_ASSERT(me_inArgs.supports(MEB::IN_ARG_alpha));
137 TEUCHOS_ASSERT(me_inArgs.supports(MEB::IN_ARG_beta));
138 }
139 TEUCHOS_ASSERT(me_outArgs.supports(MEB::OUT_ARG_DfDp, p_index_).supports(MEB::DERIV_MV_JACOBIAN_FORM));
140}
141
142template <typename Scalar>
143Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
144CombinedForwardSensitivityModelEvaluator<Scalar>::
145get_p_space(int p) const
146{
147 return model_->get_p_space(p);
148}
149
150template <typename Scalar>
151Teuchos::RCP<const Teuchos::Array<std::string> >
152CombinedForwardSensitivityModelEvaluator<Scalar>::
153get_p_names(int p) const
154{
155 return model_->get_p_names(p);
156}
157
158template <typename Scalar>
159Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
165
166template <typename Scalar>
167Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
173
174template <typename Scalar>
175Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
176CombinedForwardSensitivityModelEvaluator<Scalar>::
177get_g_space(int j) const
178{
179 if (g_index_ >= 0) {
180 if (j == 0)
181 return dgdp_space_;
182 else if (j == 1)
183 return model_->get_g_space(g_index_);
184 }
185 return model_->get_g_space(j-g_offset_);
186}
187
188template <typename Scalar>
189Teuchos::ArrayView<const std::string>
190CombinedForwardSensitivityModelEvaluator<Scalar>::
191get_g_names(int j) const
192{
193 if (g_index_ >= 0) {
194 if (j == 0) {
195 Teuchos::Array<std::string> names = model_->get_g_names(g_index_);
196 for (int i=0; i<names.size(); ++i)
197 names[i] = names[i] + "_reduced sensitivity";
198 return names();
199 }
200 else if (j == 1)
201 return model_->get_g_names(g_index_);
202 }
203 return model_->get_g_names(j-g_offset_);
204}
205
206template <typename Scalar>
207Teuchos::RCP<Thyra::LinearOpBase<Scalar> >
208CombinedForwardSensitivityModelEvaluator<Scalar>::
209create_W_op() const
210{
211 Teuchos::RCP<Thyra::LinearOpBase<Scalar> > op =
212 sens_solve_model_->create_W_op();
213 return Thyra::nonconstMultiVectorLinearOp(op, num_param_+1);
214}
215
216template <typename Scalar>
217Teuchos::RCP<Thyra::LinearOpBase<Scalar> >
218CombinedForwardSensitivityModelEvaluator<Scalar>::
219create_DgDx_dot_op(int j) const
220{
221 return model_->create_DgDx_dot_op(j-g_offset_);
222}
223
224template <typename Scalar>
225Teuchos::RCP<Thyra::LinearOpBase<Scalar> >
226CombinedForwardSensitivityModelEvaluator<Scalar>::
227create_DgDx_op(int j) const
228{
229 return model_->create_DgDx_op(j-g_offset_);
230}
231
232template <typename Scalar>
233Teuchos::RCP<Thyra::LinearOpBase<Scalar> >
234CombinedForwardSensitivityModelEvaluator<Scalar>::
235create_DgDp_op(int j, int l) const
236{
237 return model_->create_DgDp_op(j-g_offset_,l);
238}
239
240template <typename Scalar>
241Teuchos::RCP<const Thyra::LinearOpWithSolveFactoryBase<Scalar> >
242CombinedForwardSensitivityModelEvaluator<Scalar>::
243get_W_factory() const
244{
245 Teuchos::RCP<const Thyra::LinearOpWithSolveFactoryBase<Scalar> > factory =
246 sens_solve_model_->get_W_factory();
247 if (factory == Teuchos::null)
248 return Teuchos::null; // model_ doesn't support W_factory
249
250 return Thyra::multiVectorLinearOpWithSolveFactory(factory,
251 f_dfdp_space_,
252 x_dxdp_space_);
253}
254
255template <typename Scalar>
256Thyra::ModelEvaluatorBase::InArgs<Scalar>
262
263template <typename Scalar>
264Thyra::ModelEvaluatorBase::InArgs<Scalar>
265CombinedForwardSensitivityModelEvaluator<Scalar>::
266getNominalValues() const
267{
268 typedef Thyra::ModelEvaluatorBase MEB;
269 typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;
270 using Teuchos::RCP;
271 using Teuchos::rcp_dynamic_cast;
272 using Teuchos::Range1D;
273
274 MEB::InArgs<Scalar> me_nominal = model_->getNominalValues();
275 MEB::InArgs<Scalar> nominal = this->createInArgs();
276
277 const Teuchos::Range1D rng(1,num_param_);
278 const Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
279
280 // Set initial x. If dxdp_init == null, set initial dx/dp = 0
281 RCP< const Thyra::VectorBase<Scalar> > me_x = me_nominal.get_x();
282 if (me_x != Teuchos::null) {
283 RCP< Thyra::VectorBase<Scalar> > x = Thyra::createMember(*x_dxdp_space_);
284 RCP<DMVPV> x_dxdp = rcp_dynamic_cast<DMVPV>(x,true);
285 Thyra::assign(x_dxdp->getNonconstMultiVector()->col(0).ptr(), *me_x);
286 if (dxdp_init_ == Teuchos::null)
287 Thyra::assign(x_dxdp->getNonconstMultiVector()->subView(rng).ptr(),
288 zero);
289 else
290 Thyra::assign(x_dxdp->getNonconstMultiVector()->subView(rng).ptr(),
291 *dxdp_init_);
292 nominal.set_x(x);
293 }
294
295 // Set initial xdot. If dx_dotdp_init == null, set initial dxdot/dp = 0
296 RCP< const Thyra::VectorBase<Scalar> > me_xdot;
297 if (me_nominal.supports(MEB::IN_ARG_x_dot))
298 me_xdot = me_nominal.get_x_dot();
299 if (me_xdot != Teuchos::null) {
300 RCP< Thyra::VectorBase<Scalar> > xdot = Thyra::createMember(*x_dxdp_space_);
301 RCP<DMVPV> xdot_dxdp = rcp_dynamic_cast<DMVPV>(xdot,true);
302 Thyra::assign(xdot_dxdp->getNonconstMultiVector()->col(0).ptr(), *me_xdot);
303 if (dx_dotdp_init_ == Teuchos::null)
304 Thyra::assign(xdot_dxdp->getNonconstMultiVector()->subView(rng).ptr(),
305 zero);
306 else
307 Thyra::assign(xdot_dxdp->getNonconstMultiVector()->subView(rng).ptr(),
308 *dx_dotdp_init_);
309 nominal.set_x_dot(xdot);
310 }
311
312 // Set initial xdotdot. If dx_dotdotdp_init == null, set initial dxdotdot/dp = 0
313 RCP< const Thyra::VectorBase<Scalar> > me_xdotdot;
314 if (me_nominal.supports(MEB::IN_ARG_x_dot_dot))
315 me_xdotdot = me_nominal.get_x_dot_dot();
316 if (me_xdotdot != Teuchos::null) {
317 RCP< Thyra::VectorBase<Scalar> > xdotdot =
318 Thyra::createMember(*x_dxdp_space_);
319 RCP<DMVPV> xdotdot_dxdp = rcp_dynamic_cast<DMVPV>(xdotdot,true);
320 Thyra::assign(xdotdot_dxdp->getNonconstMultiVector()->col(0).ptr(), *me_xdotdot);
321 if (dx_dotdotdp_init_ == Teuchos::null)
322 Thyra::assign(xdotdot_dxdp->getNonconstMultiVector()->subView(rng).ptr(),
323 zero);
324 else
325 Thyra::assign(xdotdot_dxdp->getNonconstMultiVector()->subView(rng).ptr(),
326 *dx_dotdotdp_init_);
327 nominal.set_x_dot_dot(xdotdot);
328 }
329
330 const int np = model_->Np();
331 for (int i=0; i<np; ++i)
332 nominal.set_p(i, me_nominal.get_p(i));
333
334 return nominal;
335}
336
337template <typename Scalar>
338Thyra::ModelEvaluatorBase::OutArgs<Scalar>
344
345template <typename Scalar>
346void
347CombinedForwardSensitivityModelEvaluator<Scalar>::
348evalModelImpl(const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs,
349 const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs) const
350{
351 typedef Thyra::ModelEvaluatorBase MEB;
352 typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;
353 using Teuchos::RCP;
354 using Teuchos::rcp_dynamic_cast;
355 using Teuchos::Range1D;
356
357 const bool use_tangent = use_dfdp_as_tangent_ || use_dgdp_as_tangent_;
358
359 // setup input arguments for model
360 RCP< const Thyra::VectorBase<Scalar> > x, xdot, xdotdot;
361 RCP< const Thyra::MultiVectorBase<Scalar> > dxdp, dxdotdp, dxdotdotdp;
362 RCP< const Thyra::VectorBase<Scalar> > dxdp_vec, dxdotdp_vec, dxdotdotdp_vec;
363 MEB::InArgs<Scalar> me_inArgs = model_->getNominalValues();
364 RCP<const DMVPV> x_dxdp = rcp_dynamic_cast<const DMVPV>(inArgs.get_x(),true);
365 x = x_dxdp->getMultiVector()->col(0);
366 dxdp = x_dxdp->getMultiVector()->subView(Range1D(1,num_param_));
367 me_inArgs.set_x(x);
368 if (use_tangent) {
369 dxdp_vec = Thyra::multiVectorProductVector(dxdp_space_, dxdp);
370 if (use_dfdp_as_tangent_)
371 me_inArgs.set_p(x_tangent_index_, dxdp_vec);
372 }
373 if (me_inArgs.supports(MEB::IN_ARG_x_dot)) {
374 if (inArgs.get_x_dot() != Teuchos::null) {
375 RCP<const DMVPV> xdot_dxdotdp =
376 rcp_dynamic_cast<const DMVPV>(inArgs.get_x_dot(),true);
377 xdot = xdot_dxdotdp->getMultiVector()->col(0);
378 dxdotdp = xdot_dxdotdp->getMultiVector()->subView(Range1D(1,num_param_));
379 me_inArgs.set_x_dot(xdot);
380 if (use_tangent) {
381 dxdotdp_vec = Thyra::multiVectorProductVector(dxdp_space_, dxdotdp);
382 if (use_dfdp_as_tangent_)
383 me_inArgs.set_p(xdot_tangent_index_, dxdotdp_vec);
384 }
385 }
386 else // clear out xdot if it was set in nominalValues
387 me_inArgs.set_x_dot(Teuchos::null);
388 }
389 if (me_inArgs.supports(MEB::IN_ARG_x_dot_dot)) {
390 if (inArgs.get_x_dot_dot() != Teuchos::null) {
391 RCP<const DMVPV> xdotdot_dxdotdotdp =
392 rcp_dynamic_cast<const DMVPV>(inArgs.get_x_dot_dot(),true);
393 xdotdot = xdotdot_dxdotdotdp->getMultiVector()->col(0);
394 dxdotdotdp = xdotdot_dxdotdotdp->getMultiVector()->subView(Range1D(1,num_param_));
395 me_inArgs.set_x_dot_dot(xdotdot);
396 if (use_tangent) {
397 dxdotdotdp_vec =
398 Thyra::multiVectorProductVector(dxdp_space_, dxdotdotdp);
399 if (use_dfdp_as_tangent_)
400 me_inArgs.set_p(xdotdot_tangent_index_, dxdotdotdp_vec);
401 }
402 }
403 else // clear out xdotdot if it was set in nominalValues
404 me_inArgs.set_x_dot_dot(Teuchos::null);
405 }
406 if (me_inArgs.supports(MEB::IN_ARG_t))
407 me_inArgs.set_t(inArgs.get_t());
408 if (me_inArgs.supports(MEB::IN_ARG_alpha))
409 me_inArgs.set_alpha(inArgs.get_alpha());
410 if (me_inArgs.supports(MEB::IN_ARG_beta))
411 me_inArgs.set_beta(inArgs.get_beta());
412 if (me_inArgs.supports(MEB::IN_ARG_W_x_dot_dot_coeff))
413 me_inArgs.set_W_x_dot_dot_coeff(inArgs.get_W_x_dot_dot_coeff());
414
415 // Set parameters -- be careful to only set ones that were set in our
416 // inArgs to not null out any specified through nominalValues or
417 // dx/dp above
418 const int np = me_inArgs.Np();
419 for (int i=0; i<np; ++i) {
420 if (inArgs.get_p(i) != Teuchos::null)
421 if (!use_tangent ||
422 (use_tangent && i != x_tangent_index_ &&
423 i != xdot_tangent_index_ &&
424 i != xdotdot_tangent_index_ ))
425 me_inArgs.set_p(i, inArgs.get_p(i));
426 }
427
428 // setup output arguments for model
429 RCP< Thyra::VectorBase<Scalar> > f;
430 RCP< Thyra::MultiVectorBase<Scalar> > dfdp;
431 MEB::OutArgs<Scalar> me_outArgs = model_->createOutArgs();
432 if (outArgs.get_f() != Teuchos::null) {
433 RCP<DMVPV> f_dfdp = rcp_dynamic_cast<DMVPV>(outArgs.get_f(),true);
434 f = f_dfdp->getNonconstMultiVector()->col(0);
435 dfdp = f_dfdp->getNonconstMultiVector()->subView(Range1D(1,num_param_));
436 me_outArgs.set_f(f);
437 me_outArgs.set_DfDp(p_index_, dfdp);
438 }
439 if (outArgs.supports(MEB::OUT_ARG_W_op) &&
440 outArgs.get_W_op() != Teuchos::null &&
441 model_.ptr() == sens_solve_model_.ptr()) {
442 RCP< Thyra::LinearOpBase<Scalar> > op = outArgs.get_W_op();
443 RCP< Thyra::MultiVectorLinearOp<Scalar> > mv_op =
444 rcp_dynamic_cast< Thyra::MultiVectorLinearOp<Scalar> >(op,true);
445 me_outArgs.set_W_op(mv_op->getNonconstLinearOp());
446 }
447 for (int j=g_offset_; j<outArgs.Ng(); ++j) {
448 if (!me_outArgs.supports(MEB::OUT_ARG_DgDx_dot,j-g_offset_).none())
449 me_outArgs.set_DgDx_dot(j-g_offset_, outArgs.get_DgDx_dot(j));
450 if (!me_outArgs.supports(MEB::OUT_ARG_DgDx,j-g_offset_).none())
451 me_outArgs.set_DgDx(j-g_offset_, outArgs.get_DgDx(j));
452 for (int l=0; l<outArgs.Np(); ++l)
453 if (!me_outArgs.supports(MEB::OUT_ARG_DgDp,j-g_offset_,l).none())
454 me_outArgs.set_DgDp(j-g_offset_, l, outArgs.get_DgDp(j,l));
455 }
456 if (g_index_ >= 0 && outArgs.get_g(1) != Teuchos::null)
457 me_outArgs.set_g(g_index_, outArgs.get_g(1));
458
459 // build residual and jacobian
460 model_->evalModel(me_inArgs, me_outArgs);
461
462 // Compute W_op separately if we have a separate sensitivity solve ME
463 if (outArgs.supports(MEB::OUT_ARG_W_op) &&
464 outArgs.get_W_op() != Teuchos::null &&
465 model_.ptr() != sens_solve_model_.ptr()) {
466 MEB::OutArgs<Scalar> sens_me_outArgs = sens_solve_model_->createOutArgs();
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 sens_me_outArgs.set_W_op(mv_op->getNonconstLinearOp());
471 sens_solve_model_->evalModel(me_inArgs, sens_me_outArgs);
472 }
473
474 // Compute (df/dx) * (dx/dp) + (df/dxdot) * (dxdot/dp) + (df/dxdotdot) * (dxdotdot/dp) + (df/dp)
475 // if the underlying ME doesn't already do this. This requires computing
476 // df/dx, df/dxdot, df/dxdotdot as separate operators
477 if (!use_dfdp_as_tangent_) {
478 if (dxdp != Teuchos::null && dfdp != Teuchos::null) {
479 if (my_dfdx_ == Teuchos::null)
480 my_dfdx_ = sens_residual_model_->create_W_op();
481 MEB::OutArgs<Scalar> meo = sens_residual_model_->createOutArgs();
482 meo.set_W_op(my_dfdx_);
483 if (me_inArgs.supports(MEB::IN_ARG_alpha))
484 me_inArgs.set_alpha(0.0);
485 if (me_inArgs.supports(MEB::IN_ARG_beta))
486 me_inArgs.set_beta(1.0);
487 if (me_inArgs.supports(MEB::IN_ARG_W_x_dot_dot_coeff))
488 me_inArgs.set_W_x_dot_dot_coeff(0.0);
489 sens_residual_model_->evalModel(me_inArgs, meo);
490 my_dfdx_->apply(Thyra::NOTRANS, *dxdp, dfdp.ptr(), Scalar(1.0), Scalar(1.0));
491 }
492 if (dxdotdp != Teuchos::null && dfdp != Teuchos::null) {
493 if (my_dfdxdot_ == Teuchos::null)
494 my_dfdxdot_ = sens_residual_model_->create_W_op();
495 MEB::OutArgs<Scalar> meo = sens_residual_model_->createOutArgs();
496 meo.set_W_op(my_dfdxdot_);
497 if (me_inArgs.supports(MEB::IN_ARG_alpha))
498 me_inArgs.set_alpha(1.0);
499 if (me_inArgs.supports(MEB::IN_ARG_beta))
500 me_inArgs.set_beta(0.0);
501 if (me_inArgs.supports(MEB::IN_ARG_W_x_dot_dot_coeff))
502 me_inArgs.set_W_x_dot_dot_coeff(0.0);
503 sens_residual_model_->evalModel(me_inArgs, meo);
504 my_dfdxdot_->apply(Thyra::NOTRANS, *dxdotdp, dfdp.ptr(), Scalar(1.0), Scalar(1.0));
505 }
506 if (dxdotdotdp != Teuchos::null && dfdp != Teuchos::null) {
507 if (my_dfdxdotdot_ == Teuchos::null)
508 my_dfdxdotdot_ = sens_residual_model_->create_W_op();
509 MEB::OutArgs<Scalar> meo = sens_residual_model_->createOutArgs();
510 meo.set_W_op(my_dfdxdotdot_);
511 if (me_inArgs.supports(MEB::IN_ARG_alpha))
512 me_inArgs.set_alpha(0.0);
513 if (me_inArgs.supports(MEB::IN_ARG_beta))
514 me_inArgs.set_beta(0.0);
515 if (me_inArgs.supports(MEB::IN_ARG_W_x_dot_dot_coeff))
516 me_inArgs.set_W_x_dot_dot_coeff(1.0);
517 sens_residual_model_->evalModel(me_inArgs, meo);
518 my_dfdxdotdot_->apply(Thyra::NOTRANS, *dxdotdotdp, dfdp.ptr(), Scalar(1.0), Scalar(1.0));
519 }
520 }
521
522 // Reduced dg/dp = dg/dx*dx/dp + dg/dp
523 if (g_index_ >= 0 && outArgs.get_g(0) != Teuchos::null) {
524 MEB::OutArgs<Scalar> meo = model_->createOutArgs();
525 RCP<Thyra::MultiVectorBase<Scalar> > dgdp =
526 rcp_dynamic_cast<DMVPV>(outArgs.get_g(0),true)->getNonconstMultiVector();
527 RCP<Thyra::MultiVectorBase<Scalar> > dgdp_trans;
528 MEB::DerivativeSupport dgdp_support =
529 meo.supports(MEB::OUT_ARG_DgDp, g_index_, p_index_);
530 if (dgdp_support.supports(MEB::DERIV_MV_JACOBIAN_FORM))
531 meo.set_DgDp(g_index_, p_index_,
532 MEB::Derivative<Scalar>(dgdp, MEB::DERIV_MV_JACOBIAN_FORM));
533 else if (dgdp_support.supports(MEB::DERIV_MV_GRADIENT_FORM)) {
534 dgdp_trans = createMembers(model_->get_p_space(p_index_),
535 num_response_);
536 meo.set_DgDp(g_index_, p_index_,
537 MEB::Derivative<Scalar>(dgdp_trans, MEB::DERIV_MV_GRADIENT_FORM));
538 }
539 else
540 TEUCHOS_TEST_FOR_EXCEPTION(
541 true, std::logic_error,
542 "Operator form of dg/dp not supported for reduced sensitivity");
543
544 if (!use_dgdp_as_tangent_ || dxdp == Teuchos::null) {
545 MEB::DerivativeSupport dgdx_support =
546 meo.supports(MEB::OUT_ARG_DgDx, g_index_);
547 if (dgdx_support.supports(MEB::DERIV_LINEAR_OP)) {
548 if (my_dgdx_ == Teuchos::null)
549 my_dgdx_ = model_->create_DgDx_op(g_index_);
550 meo.set_DgDx(g_index_, my_dgdx_);
551 }
552 else if (dgdx_support.supports(MEB::DERIV_MV_GRADIENT_FORM)) {
553 if (my_dgdx_mv_ == Teuchos::null)
554 my_dgdx_mv_ = createMembers(model_->get_g_space(g_index_), num_param_);
555 meo.set_DgDx(g_index_,
556 MEB::Derivative<Scalar>(my_dgdx_mv_,
557 MEB::DERIV_MV_GRADIENT_FORM));
558 }
559 else
560 TEUCHOS_TEST_FOR_EXCEPTION(
561 true, std::logic_error,
562 "Jacobian form of dg/dx not supported for reduced sensitivity");
563 }
564
565 // Clear dx/dp, dxdot/dp, dxdotdot/dp from inArgs if set and we are not
566 // using dg/dp as the tangent
567 if (!use_dgdp_as_tangent_ && use_dfdp_as_tangent_ &&
568 dxdp_vec != Teuchos::null)
569 me_inArgs.set_p(x_tangent_index_, Teuchos::null);
570 if (!use_dgdp_as_tangent_ && use_dfdp_as_tangent_ &&
571 dxdotdp_vec != Teuchos::null)
572 me_inArgs.set_p(xdot_tangent_index_, Teuchos::null);
573 if (!use_dgdp_as_tangent_ && use_dfdp_as_tangent_ &&
574 dxdotdotdp_vec != Teuchos::null)
575 me_inArgs.set_p(xdotdot_tangent_index_, Teuchos::null);
576
577 // Set dx/dp, dxdot/dp, dxdodot/dp if necessary
578 if (use_dgdp_as_tangent_ && !use_dfdp_as_tangent_ &&
579 dxdp_vec != Teuchos::null)
580 me_inArgs.set_p(x_tangent_index_, dxdp_vec);
581 if (use_dgdp_as_tangent_ && !use_dfdp_as_tangent_ &&
582 dxdotdp_vec != Teuchos::null)
583 me_inArgs.set_p(xdot_tangent_index_, dxdotdp_vec);
584 if (use_dgdp_as_tangent_ && !use_dfdp_as_tangent_ &&
585 dxdotdotdp_vec != Teuchos::null)
586 me_inArgs.set_p(xdotdot_tangent_index_, dxdotdotdp_vec);
587
588 model_->evalModel(me_inArgs, meo);
589
590 // transpose reduced dg/dp if necessary
591 if (dgdp_trans != Teuchos::null) {
592 Thyra::DetachedMultiVectorView<Scalar> dgdp_view(*dgdp);
593 Thyra::DetachedMultiVectorView<Scalar> dgdp_trans_view(*dgdp_trans);
594 for (int i=0; i<num_param_; ++i)
595 for (int j=0; j<num_response_; ++j)
596 dgdp_view(j,i) = dgdp_trans_view(i,j);
597 }
598
599 // Compute (dg/dx) * (dx/dp) + (dg/dp) if the underlying ME doesn't already
600 // do this.
601 if (!use_dgdp_as_tangent_ && dxdp != Teuchos::null) {
602 MEB::DerivativeSupport dgdx_support =
603 meo.supports(MEB::OUT_ARG_DgDx, g_index_);
604 if (dgdx_support.supports(MEB::DERIV_LINEAR_OP))
605 my_dgdx_->apply(Thyra::NOTRANS, *dxdp, dgdp.ptr(), Scalar(1.0), Scalar(1.0));
606 else if (dgdx_support.supports(MEB::DERIV_MV_GRADIENT_FORM))
607 my_dgdx_mv_->apply(Thyra::TRANS, *dxdp, dgdp.ptr(), Scalar(1.0), Scalar(1.0));
608 else
609 TEUCHOS_TEST_FOR_EXCEPTION(
610 true, std::logic_error,
611 "Jacobian form of dg/dx not supported for reduced sensitivity");
612 }
613 }
614}
615
616template<class Scalar>
617Teuchos::RCP<const Teuchos::ParameterList>
618CombinedForwardSensitivityModelEvaluator<Scalar>::
619getValidParameters()
620{
621 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
622 pl->set<bool>("Use DfDp as Tangent", false);
623 pl->set<bool>("Use DgDp as Tangent", false);
624 pl->set<int>("Sensitivity Parameter Index", 0);
625 pl->set<int>("Response Function Index", -1);
626 pl->set<int>("Sensitivity X Tangent Index", 1);
627 pl->set<int>("Sensitivity X-Dot Tangent Index", 2);
628 pl->set<int>("Sensitivity X-Dot-Dot Tangent Index", 3);
629 return pl;
630}
631
632} // namespace Tempus
633
634#endif
Thyra_MultiVectorLinearOpWithSolveFactory.hpp
Thyra_MultiVectorLinearOp.hpp
Thyra_ReuseLinearOpWithSolveFactory.hpp
Tempus::CombinedForwardSensitivityModelEvaluator::get_p_space
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_p_space(int p) const
Definition Tempus_CombinedForwardSensitivityModelEvaluator_impl.hpp:145
Tempus::CombinedForwardSensitivityModelEvaluator::get_p_names
Teuchos::RCP< const Teuchos::Array< std::string > > get_p_names(int p) const
Definition Tempus_CombinedForwardSensitivityModelEvaluator_impl.hpp:153
Tempus::CombinedForwardSensitivityModelEvaluator::get_g_space
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_g_space(int j) const
Definition Tempus_CombinedForwardSensitivityModelEvaluator_impl.hpp:177
Tempus::CombinedForwardSensitivityModelEvaluator::create_DgDx_op
Teuchos::RCP< Thyra::LinearOpBase< Scalar > > create_DgDx_op(int j) const
Definition Tempus_CombinedForwardSensitivityModelEvaluator_impl.hpp:227
Tempus::CombinedForwardSensitivityModelEvaluator::get_W_factory
Teuchos::RCP< const Thyra::LinearOpWithSolveFactoryBase< Scalar > > get_W_factory() const
Definition Tempus_CombinedForwardSensitivityModelEvaluator_impl.hpp:243
Tempus::CombinedForwardSensitivityModelEvaluator::evalModelImpl
void evalModelImpl(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs, const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &outArgs) const
Definition Tempus_CombinedForwardSensitivityModelEvaluator_impl.hpp:348
Tempus::CombinedForwardSensitivityModelEvaluator::create_DgDp_op
Teuchos::RCP< Thyra::LinearOpBase< Scalar > > create_DgDp_op(int j, int l) const
Definition Tempus_CombinedForwardSensitivityModelEvaluator_impl.hpp:235
Tempus::CombinedForwardSensitivityModelEvaluator::create_DgDx_dot_op
Teuchos::RCP< Thyra::LinearOpBase< Scalar > > create_DgDx_dot_op(int j) const
Definition Tempus_CombinedForwardSensitivityModelEvaluator_impl.hpp:219
Tempus::CombinedForwardSensitivityModelEvaluator::get_f_space
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_f_space() const
Definition Tempus_CombinedForwardSensitivityModelEvaluator_impl.hpp:169
Tempus::CombinedForwardSensitivityModelEvaluator::CombinedForwardSensitivityModelEvaluator
CombinedForwardSensitivityModelEvaluator(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 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.
Definition Tempus_CombinedForwardSensitivityModelEvaluator_impl.hpp:23
Tempus::CombinedForwardSensitivityModelEvaluator::get_x_space
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_x_space() const
Definition Tempus_CombinedForwardSensitivityModelEvaluator_impl.hpp:161
Thyra::MultiVectorLinearOpWithSolveFactory::multiVectorLinearOpWithSolveFactory
RCP< MultiVectorLinearOpWithSolveFactory< Scalar > > multiVectorLinearOpWithSolveFactory(const RCP< const LinearOpWithSolveFactoryBase< Scalar > > &lowsf, const RCP< const DefaultMultiVectorProductVectorSpace< Scalar > > &multiVecRange, const RCP< const DefaultMultiVectorProductVectorSpace< Scalar > > &multiVecDomain)
Nonmember constructor.
Definition Thyra_MultiVectorLinearOpWithSolveFactory.hpp:226
Thyra::MultiVectorLinearOp::nonconstMultiVectorLinearOp
RCP< MultiVectorLinearOp< Scalar > > nonconstMultiVectorLinearOp(const RCP< LinearOpBase< Scalar > > &op, const RCP< const DefaultMultiVectorProductVectorSpace< Scalar > > &multiVecRange, const RCP< const DefaultMultiVectorProductVectorSpace< Scalar > > &multiVecDomain)
Nonmember constructor function.
Definition Thyra_MultiVectorLinearOp.hpp:275
Generated by doxygen 1.13.2