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