Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
Tempus_AdjointAuxSensitivityModelEvaluator_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_AdjointAuxSensitivityModelEvaluator_impl_hpp
10#define Tempus_AdjointAuxSensitivityModelEvaluator_impl_hpp
11
16#include "Thyra_DefaultBlockedLinearOp.hpp"
18#include "Thyra_VectorStdOps.hpp"
19#include "Thyra_MultiVectorStdOps.hpp"
20
21namespace Tempus {
22
23template <typename Scalar>
26 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> > & model,
27 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> > & adjoint_model,
28 const Scalar& t_init,
29 const Scalar& t_final,
30 const Teuchos::RCP<const Teuchos::ParameterList>& pList) :
31 model_(model),
32 adjoint_model_(adjoint_model),
33 t_init_(t_init),
34 t_final_(t_final),
35 mass_matrix_is_computed_(false),
36 t_interp_(Teuchos::ScalarTraits<Scalar>::rmax())
37{
38 using Teuchos::RCP;
39 using Teuchos::Array;
40 using Thyra::VectorSpaceBase;
41 typedef Thyra::ModelEvaluatorBase MEB;
42
43 // Set parameters
44 Teuchos::RCP<Teuchos::ParameterList> pl =
45 Teuchos::rcp(new Teuchos::ParameterList);
46 if (pList != Teuchos::null)
47 *pl = *pList;
48 pl->validateParametersAndSetDefaults(*this->getValidParameters());
49 mass_matrix_is_constant_ = pl->get<bool>("Mass Matrix Is Constant");
50 mass_matrix_is_identity_ = pl->get<bool>("Mass Matrix Is Identity");
51 p_index_ = pl->get<int>("Sensitivity Parameter Index", 0);
52 g_index_ = pl->get<int>("Response Function Index", 0);
53 num_adjoint_ = model_->get_g_space(g_index_)->dim();
54
55 // We currently do not support a non-constant mass matrix
56 TEUCHOS_TEST_FOR_EXCEPTION(
57 mass_matrix_is_constant_ == false, std::logic_error,
58 "AdjointAuxSensitivityModelEvaluator currently does not support " <<
59 "non-constant mass matrix df/dx_dot!");
60
62 Thyra::multiVectorProductVectorSpace(model_->get_f_space(), num_adjoint_);
64 Thyra::multiVectorProductVectorSpace(model_->get_x_space(), num_adjoint_);
66 Thyra::multiVectorProductVectorSpace(model_->get_p_space(p_index_),
68 Array< RCP<const VectorSpaceBase<Scalar> > > x_spaces(2);
69 Array< RCP<const VectorSpaceBase<Scalar> > > f_spaces(2);
70 x_spaces[0] = adjoint_space_;
71 x_spaces[1] = response_space_;
72 f_spaces[0] = residual_space_;
73 f_spaces[1] = response_space_;
74 x_prod_space_ = Thyra::productVectorSpace(x_spaces());
75 f_prod_space_ = Thyra::productVectorSpace(f_spaces());
76
77 // forward and adjoint models must support same InArgs
78 MEB::InArgs<Scalar> me_inArgs = model_->createInArgs();
79 me_inArgs.assertSameSupport(adjoint_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 inArgs.setSupports(MEB::IN_ARG_alpha);
88 inArgs.setSupports(MEB::IN_ARG_beta);
89
90 // Support additional parameters for x and xdot
91 inArgs.set_Np(me_inArgs.Np());
92 prototypeInArgs_ = inArgs;
93
94 MEB::OutArgs<Scalar> me_outArgs = model_->createOutArgs();
95 MEB::OutArgs<Scalar> adj_me_outArgs = adjoint_model_->createOutArgs();
96 MEB::OutArgsSetup<Scalar> outArgs;
97 outArgs.setModelEvalDescription(this->description());
98 outArgs.set_Np_Ng(me_inArgs.Np(),0);
99 outArgs.setSupports(MEB::OUT_ARG_f);
100 outArgs.setSupports(MEB::OUT_ARG_W_op);
101 prototypeOutArgs_ = outArgs;
102
103 // Adjoint ME must support W_op to define adjoint ODE/DAE.
104 // Must support alpha, beta if it suports x_dot
105 TEUCHOS_ASSERT(me_inArgs.supports(MEB::IN_ARG_x));
106 TEUCHOS_ASSERT(adj_me_outArgs.supports(MEB::OUT_ARG_W_op));
107 if (me_inArgs.supports(MEB::IN_ARG_x_dot)) {
108 TEUCHOS_ASSERT(me_inArgs.supports(MEB::IN_ARG_alpha));
109 TEUCHOS_ASSERT(me_inArgs.supports(MEB::IN_ARG_beta));
110 }
111}
112
113template <typename Scalar>
114void
116setFinalTime(const Scalar t_final)
117{
118 t_final_ = t_final;
119}
120
121template <typename Scalar>
122void
125 const Teuchos::RCP<const Tempus::SolutionHistory<Scalar> >& sh)
126{
127 sh_ = sh;
128 t_interp_ = Teuchos::ScalarTraits<Scalar>::rmax();
129 forward_state_ = Teuchos::null;
130}
131
132template <typename Scalar>
133Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
135get_p_space(int p) const
136{
137 TEUCHOS_ASSERT(p < model_->Np());
138 return model_->get_p_space(p);
139}
140
141template <typename Scalar>
142Teuchos::RCP<const Teuchos::Array<std::string> >
144get_p_names(int p) const
145{
146 TEUCHOS_ASSERT(p < model_->Np());
147 return model_->get_p_names(p);
148}
149
150template <typename Scalar>
151Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
153get_x_space() const
154{
155 return x_prod_space_;
156}
157
158template <typename Scalar>
159Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
161get_f_space() const
162{
163 return f_prod_space_;
164}
165
166template <typename Scalar>
167Teuchos::RCP<Thyra::LinearOpBase<Scalar> >
169create_W_op() const
170{
171 using Teuchos::RCP;
173
174 RCP<LinearOpBase<Scalar> > adjoint_op = adjoint_model_->create_W_op();
175 RCP<LinearOpBase<Scalar> > mv_adjoint_op =
176 Thyra::nonconstMultiVectorLinearOp(adjoint_op, num_adjoint_);
177 RCP<const Thyra::VectorSpaceBase<Scalar> > g_space = response_space_;
178 RCP<LinearOpBase<Scalar> > g_op =
179 Thyra::scaledIdentity(g_space, Scalar(1.0));
180 RCP<LinearOpBase<Scalar> > null_op;
181 return nonconstBlock2x2(mv_adjoint_op, null_op, null_op, g_op);
182}
183
184template <typename Scalar>
185Teuchos::RCP<const Thyra::LinearOpWithSolveFactoryBase<Scalar> >
187get_W_factory() const
188{
189 using Teuchos::RCP;
190 using Teuchos::rcp_dynamic_cast;
192
193 RCP<const LOWSFB > alowsfb = adjoint_model_->get_W_factory();
194 if (alowsfb == Teuchos::null)
195 return Teuchos::null; // model_ doesn't support W_factory
196
197 RCP<const LOWSFB > mv_alowsfb =
198 Thyra::multiVectorLinearOpWithSolveFactory(alowsfb, residual_space_,
199 adjoint_space_);
200
201 RCP<const Thyra::VectorSpaceBase<Scalar> > g_space = response_space_;
202 RCP<const LOWSFB > g_lowsfb =
203 Thyra::scaledIdentitySolveFactory(g_space, Scalar(1.0));
204
205 Teuchos::Array< RCP<const LOWSFB > > lowsfbs(2);
206 lowsfbs[0] = mv_alowsfb;
207 lowsfbs[1] = g_lowsfb;
208 return Thyra::blockedTriangularLinearOpWithSolveFactory(lowsfbs);
209}
210
211template <typename Scalar>
212Thyra::ModelEvaluatorBase::InArgs<Scalar>
214createInArgs() const
215{
216 return prototypeInArgs_;
217}
218
219template <typename Scalar>
220Thyra::ModelEvaluatorBase::InArgs<Scalar>
222getNominalValues() const
223{
224 typedef Thyra::ModelEvaluatorBase MEB;
225 using Teuchos::RCP;
226 using Teuchos::rcp_dynamic_cast;
227
228 MEB::InArgs<Scalar> me_nominal = model_->getNominalValues();
229 MEB::InArgs<Scalar> nominal = this->createInArgs();
230
231 const Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
232
233 // Set initial x, x_dot
234 RCP< Thyra::VectorBase<Scalar> > x = Thyra::createMember(*x_prod_space_);
235 Thyra::assign(x.ptr(), zero);
236 nominal.set_x(x);
237
238 if (me_nominal.supports(MEB::IN_ARG_x_dot)) {
239 RCP< Thyra::VectorBase<Scalar> > x_dot =
240 Thyra::createMember(*x_prod_space_);
241 Thyra::assign(x_dot.ptr(), zero);
242 nominal.set_x_dot(x_dot);
243 }
244
245 const int np = model_->Np();
246 for (int i=0; i<np; ++i)
247 nominal.set_p(i, me_nominal.get_p(i));
248
249 return nominal;
250}
251
252template <typename Scalar>
253Thyra::ModelEvaluatorBase::OutArgs<Scalar>
255createOutArgsImpl() const
256{
257 return prototypeOutArgs_;
258}
259
260template <typename Scalar>
261void
263evalModelImpl(const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs,
264 const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs) const
265{
266 typedef Thyra::ModelEvaluatorBase MEB;
267 using Teuchos::RCP;
268 using Teuchos::rcp_dynamic_cast;
269
270 // Note: adjoint_model computes the transposed W (either explicitly or
271 // implicitly. Thus we need to always call adjoint_model->evalModel()
272 // whenever computing the adjoint operator, and subsequent calls to apply()
273 // do not transpose it.
274
275 // Interpolate forward solution at supplied time, reusing previous
276 // interpolation if possible
277 TEUCHOS_ASSERT(sh_ != Teuchos::null);
278 const Scalar t = inArgs.get_t();
279 const Scalar forward_t = t_final_ + t_init_ - t;
280 if (forward_state_ == Teuchos::null || t_interp_ != t) {
281 if (forward_state_ == Teuchos::null)
282 forward_state_ = sh_->interpolateState(forward_t);
283 else
284 sh_->interpolateState(forward_t, forward_state_.get());
285 t_interp_ = t;
286 }
287
288 // setup input arguments for model
289 MEB::InArgs<Scalar> me_inArgs = model_->getNominalValues();
290 me_inArgs.set_x(forward_state_->getX());
291 if (me_inArgs.supports(MEB::IN_ARG_x_dot) &&
292 inArgs.get_x_dot() != Teuchos::null)
293 me_inArgs.set_x_dot(forward_state_->getXDot());
294 if (me_inArgs.supports(MEB::IN_ARG_t))
295 me_inArgs.set_t(forward_t);
296 const int np = me_inArgs.Np();
297 for (int i=0; i<np; ++i)
298 me_inArgs.set_p(i, inArgs.get_p(i));
299
300 // compute W
301 RCP<Thyra::LinearOpBase<Scalar> > op = outArgs.get_W_op();
302 if (op != Teuchos::null) {
303 if (me_inArgs.supports(MEB::IN_ARG_alpha))
304 me_inArgs.set_alpha(inArgs.get_alpha());
305 if (me_inArgs.supports(MEB::IN_ARG_beta))
306 me_inArgs.set_beta(inArgs.get_beta());
307
308 // Adjoint W
309 RCP<Thyra::DefaultBlockedLinearOp<Scalar> > block_op =
310 rcp_dynamic_cast<Thyra::DefaultBlockedLinearOp<Scalar> >(op,true);
311 RCP<Thyra::MultiVectorLinearOp<Scalar> > mv_adjoint_op =
312 rcp_dynamic_cast<Thyra::MultiVectorLinearOp<Scalar> >(
313 block_op->getNonconstBlock(0,0),true);
314 RCP<Thyra::LinearOpBase<Scalar> > adjoint_op =
315 mv_adjoint_op->getNonconstLinearOp();
316 MEB::OutArgs<Scalar> adj_me_outArgs = adjoint_model_->createOutArgs();
317 adj_me_outArgs.set_W_op(adjoint_op);
318 adjoint_model_->evalModel(me_inArgs, adj_me_outArgs);
319
320 // g W
321 RCP<Thyra::ScaledIdentityLinearOpWithSolve<Scalar> > si_op =
322 rcp_dynamic_cast<Thyra::ScaledIdentityLinearOpWithSolve<Scalar> >(
323 block_op->getNonconstBlock(1,1),true);
324 si_op->setScale(inArgs.get_alpha());
325 }
326
327 // Compute adjoint residual F(y):
328 // * For implicit form, F(y) = d/dt( df/dx_dot^T*y ) + df/dx^T*y
329 // * For explict form, F(y) = -df/dx^T*y
330 // For implicit form, we assume df/dx_dot is constant w.r.t. x, x_dot, and t,
331 // so the residual becomes F(y) = df/dx_dot^T*y_dot + df/dx^T*y
332 RCP<Thyra::VectorBase<Scalar> > f = outArgs.get_f();
333 if (f != Teuchos::null) {
334 RCP<const Thyra::VectorBase<Scalar> > x = inArgs.get_x().assert_not_null();
335 RCP<const DPV> prod_x = rcp_dynamic_cast<const DPV>(x,true);
336 RCP<const Thyra::VectorBase<Scalar> > adjoint_x = prod_x->getVectorBlock(0);
337 RCP<const Thyra::MultiVectorBase<Scalar> >adjoint_x_mv =
338 rcp_dynamic_cast<const DMVPV>(adjoint_x,true)->getMultiVector();
339
340 RCP<DPV> prod_f = rcp_dynamic_cast<DPV>(f,true);
341 RCP<Thyra::VectorBase<Scalar> > adjoint_f =
342 prod_f->getNonconstVectorBlock(0);
343 RCP<Thyra::MultiVectorBase<Scalar> > adjoint_f_mv =
344 rcp_dynamic_cast<DMVPV>(adjoint_f,true)->getNonconstMultiVector();
345
346 MEB::OutArgs<Scalar> adj_me_outArgs = adjoint_model_->createOutArgs();
347
348 if (my_dfdx_ == Teuchos::null)
349 my_dfdx_ = adjoint_model_->create_W_op();
350 adj_me_outArgs.set_W_op(my_dfdx_);
351 if (me_inArgs.supports(MEB::IN_ARG_alpha))
352 me_inArgs.set_alpha(0.0);
353 if (me_inArgs.supports(MEB::IN_ARG_beta))
354 me_inArgs.set_beta(1.0);
355 adjoint_model_->evalModel(me_inArgs, adj_me_outArgs);
356
357 // Explicit form residual F(y) = -df/dx^T*y
358 my_dfdx_->apply(Thyra::NOTRANS, *adjoint_x_mv, adjoint_f_mv.ptr(),
359 Scalar(-1.0), Scalar(0.0));
360
361 // Implicit form residual df/dx_dot^T*y_dot + df/dx^T*y using the second
362 // scalar argument to apply() to change the explicit term above
363 RCP<const DPV> prod_x_dot;
364 if (me_inArgs.supports(MEB::IN_ARG_x_dot)) {
365 RCP<const Thyra::VectorBase<Scalar> > x_dot = inArgs.get_x_dot();
366 if (x_dot != Teuchos::null) {
367 prod_x_dot = rcp_dynamic_cast<const DPV>(x_dot,true);
368 RCP<const Thyra::VectorBase<Scalar> > adjoint_x_dot =
369 prod_x_dot->getVectorBlock(0);
370 RCP<const Thyra::MultiVectorBase<Scalar> > adjoint_x_dot_mv =
371 rcp_dynamic_cast<const DMVPV>(adjoint_x_dot,true)->getMultiVector();
372 if (mass_matrix_is_identity_) {
373 // F = -F + y_dot
374 Thyra::V_StVpV(adjoint_f_mv.ptr(), Scalar(-1.0), *adjoint_f_mv,
375 *adjoint_x_dot_mv);
376 }
377 else {
378 if (my_dfdxdot_ == Teuchos::null)
379 my_dfdxdot_ = adjoint_model_->create_W_op();
380 if (!mass_matrix_is_constant_ || !mass_matrix_is_computed_) {
381 adj_me_outArgs.set_W_op(my_dfdxdot_);
382 me_inArgs.set_alpha(1.0);
383 me_inArgs.set_beta(0.0);
384 adjoint_model_->evalModel(me_inArgs, adj_me_outArgs);
385
386 mass_matrix_is_computed_ = true;
387 }
388 my_dfdxdot_->apply(Thyra::NOTRANS, *adjoint_x_dot_mv,
389 adjoint_f_mv.ptr(), Scalar(1.0), Scalar(-1.0));
390 }
391 }
392 }
393
394 // Compute g = z_dot - df/dp^T*y for computing the model parameter term
395 // in the adjoint sensitivity formula
396 RCP<Thyra::VectorBase<Scalar> > adjoint_g =
397 prod_f->getNonconstVectorBlock(1);
398 RCP<Thyra::MultiVectorBase<Scalar> > adjoint_g_mv =
399 rcp_dynamic_cast<DMVPV>(adjoint_g,true)->getNonconstMultiVector();
400
401 MEB::OutArgs<Scalar> me_outArgs2 = model_->createOutArgs();
402 MEB::DerivativeSupport dfdp_support =
403 me_outArgs2.supports(MEB::OUT_ARG_DfDp, p_index_);
404 Thyra::EOpTransp trans = Thyra::CONJTRANS;
405 if (dfdp_support.supports(MEB::DERIV_LINEAR_OP)) {
406 if (my_dfdp_op_ == Teuchos::null)
407 my_dfdp_op_ = model_->create_DfDp_op(p_index_);
408 me_outArgs2.set_DfDp(p_index_, MEB::Derivative<Scalar>(my_dfdp_op_));
409 trans = Thyra::CONJTRANS;
410 }
411 else if (dfdp_support.supports(MEB::DERIV_MV_JACOBIAN_FORM)) {
412 if (my_dfdp_mv_ == Teuchos::null)
413 my_dfdp_mv_ = createMembers(model_->get_f_space(),
414 model_->get_p_space(p_index_)->dim());
415 me_outArgs2.set_DfDp(p_index_,
416 MEB::Derivative<Scalar>(my_dfdp_mv_,
417 MEB::DERIV_MV_JACOBIAN_FORM));
418 my_dfdp_op_ = my_dfdp_mv_;
419 trans = Thyra::CONJTRANS;
420 }
421 else if (dfdp_support.supports(MEB::DERIV_MV_GRADIENT_FORM)) {
422 if (my_dfdp_mv_ == Teuchos::null)
423 my_dfdp_mv_ = createMembers(model_->get_p_space(p_index_),
424 model_->get_f_space()->dim());
425 me_outArgs2.set_DfDp(p_index_,
426 MEB::Derivative<Scalar>(my_dfdp_mv_,
427 MEB::DERIV_MV_GRADIENT_FORM));
428 my_dfdp_op_ = my_dfdp_mv_;
429 trans = Thyra::CONJ;
430 }
431 else
432 TEUCHOS_TEST_FOR_EXCEPTION(
433 true, std::logic_error, "Invalid df/dp support");
434 model_->evalModel(me_inArgs, me_outArgs2);
435 my_dfdp_op_->apply(trans, *adjoint_x_mv, adjoint_g_mv.ptr(),
436 Scalar(1.0), Scalar(0.0));
437
438 if (prod_x_dot != Teuchos::null) {
439 RCP<const Thyra::VectorBase<Scalar> > z_dot =
440 prod_x_dot->getVectorBlock(1);
441 RCP<const Thyra::MultiVectorBase<Scalar> > z_dot_mv =
442 rcp_dynamic_cast<const DMVPV>(z_dot,true)->getMultiVector();
443 Thyra::V_VmV(adjoint_g_mv.ptr(), *z_dot_mv, *adjoint_g_mv);
444 }
445 }
446}
447
448template<class Scalar>
449Teuchos::RCP<const Teuchos::ParameterList>
452{
453 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
454 pl->set<int>("Sensitivity Parameter Index", 0);
455 pl->set<int>("Response Function Index", 0);
456 pl->set<bool>("Mass Matrix Is Constant", true);
457 pl->set<bool>("Mass Matrix Is Identity", false);
458 return pl;
459}
460
461} // namespace Tempus
462
463#endif
void evalModelImpl(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs, const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &outArgs) const
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_p_space(int p) const
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_x_space() const
Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > adjoint_model_
Teuchos::RCP< Thyra::LinearOpBase< Scalar > > create_W_op() const
static Teuchos::RCP< const Teuchos::ParameterList > getValidParameters()
void setFinalTime(const Scalar t_final)
Set the final time from the forward evaluation.
Thyra::ModelEvaluatorBase::InArgs< Scalar > createInArgs() const
Thyra::ModelEvaluatorBase::InArgs< Scalar > getNominalValues() const
Thyra::ModelEvaluatorBase::OutArgs< Scalar > createOutArgsImpl() 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
Teuchos::RCP< const Thyra::LinearOpWithSolveFactoryBase< Scalar > > get_W_factory() const
Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > model_
Teuchos::RCP< const Teuchos::Array< std::string > > get_p_names(int p) const
AdjointAuxSensitivityModelEvaluator(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &adjoint_model, const Scalar &t_init, const Scalar &t_final, const Teuchos::RCP< const Teuchos::ParameterList > &pList=Teuchos::null)
Constructor.
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...