Rythmos - Transient Integration for Differential Equations Version of the Day
Loading...
Searching...
No Matches
Rythmos_ForwardResponseSensitivityComputer.hpp
1//@HEADER
2// ***********************************************************************
3//
4// Rythmos Package
5// Copyright (2006) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// This library is free software; you can redistribute it and/or modify
11// it under the terms of the GNU Lesser General Public License as
12// published by the Free Software Foundation; either version 2.1 of the
13// License, or (at your option) any later version.
14//
15// This library is distributed in the hope that it will be useful, but
16// WITHOUT ANY WARRANTY; without even the implied warranty of
17// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18// Lesser General Public License for more details.
19//
20// You should have received a copy of the GNU Lesser General Public
21// License along with this library; if not, write to the Free Software
22// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
23// USA
24// Questions? Contact Todd S. Coffey (tscoffe@sandia.gov)
25//
26// ***********************************************************************
27//@HEADER
28
29#ifndef RYTHMOS_FORWARD_RESPONSE_SENSITIVITY_COMPUTER_HPP
30#define RYTHMOS_FORWARD_RESPONSE_SENSITIVITY_COMPUTER_HPP
31
32
33#include "Rythmos_Types.hpp"
34#include "Thyra_ModelEvaluator.hpp"
35#include "Teuchos_VerboseObject.hpp"
36#include "Teuchos_StandardMemberCompositionMacros.hpp"
37
38
39namespace Rythmos {
40
41
48template<class Scalar>
50 : public Teuchos::VerboseObject<ForwardResponseSensitivityComputer<Scalar> >
51{
52public:
53
56
58 STANDARD_MEMBER_COMPOSITION_MEMBERS( bool, dumpSensitivities );
59
77 const RCP<const Thyra::ModelEvaluator<Scalar> > &responseFunc,
78 const Thyra::ModelEvaluatorBase::InArgs<Scalar> &basePoint,
79 const int p_index,
80 const int g_index
81 );
82
88 const RCP<const Thyra::ModelEvaluator<Scalar> > &responseFunc,
89 const Thyra::ModelEvaluatorBase::InArgs<Scalar> &basePoint
90 );
91
93 const RCP<Thyra::VectorBase<Scalar> > create_g_hat() const;
94
96 const RCP<Thyra::MultiVectorBase<Scalar> > create_D_g_hat_D_p() const;
97
110 void computeResponse(
111 const Thyra::VectorBase<Scalar> *x_dot,
112 const Thyra::VectorBase<Scalar> &x,
113 const Scalar t,
114 Thyra::VectorBase<Scalar> *g_hat
115 ) const;
116
137 const Thyra::VectorBase<Scalar> *x_dot,
138 const Thyra::MultiVectorBase<Scalar> *S_dot,
139 const Thyra::VectorBase<Scalar> &x,
140 const Thyra::MultiVectorBase<Scalar> &S,
141 const Scalar t,
142 Thyra::VectorBase<Scalar> *g_hat,
143 Thyra::MultiVectorBase<Scalar> *D_g_hat_D_p
144 ) const;
145
146private: // Data members
147
148 RCP<const Thyra::ModelEvaluator<Scalar> > responseFunc_;
149 Thyra::ModelEvaluatorBase::InArgs<Scalar> basePoint_;
150 int p_index_;
151 int g_index_;
152
153 RCP<const Thyra::VectorSpaceBase<Scalar> > p_space_;
154 RCP<const Thyra::VectorSpaceBase<Scalar> > g_space_;
155
156 bool response_func_supports_x_dot_;
157 bool response_func_supports_D_x_dot_;
158 bool response_func_supports_D_p_;
159
160 mutable Thyra::ModelEvaluatorBase::InArgs<Scalar> responseInArgs_;
161 mutable Thyra::ModelEvaluatorBase::OutArgs<Scalar> responseOutArgs_;
162
163 mutable RCP<Thyra::LinearOpBase<Scalar> > D_g_D_x_dot_;
164 mutable RCP<Thyra::LinearOpBase<Scalar> > D_g_D_x_;
165 mutable RCP<Thyra::MultiVectorBase<Scalar> > D_g_D_p_;
166
167private: // Functions
168
169 void clearCache();
170
171 void createCache(const bool computeSens) const;
172
173 void computeResponseAndSensitivityImpl(
174 const Thyra::VectorBase<Scalar> *x_dot,
175 const Thyra::MultiVectorBase<Scalar> *S_dot,
176 const Thyra::VectorBase<Scalar> &x,
177 const Thyra::MultiVectorBase<Scalar> *S,
178 const Scalar t,
179 Thyra::VectorBase<Scalar> *g_hat,
180 Thyra::MultiVectorBase<Scalar> *D_g_hat_D_p
181 ) const;
182
183};
184
185
186//
187// Implementations
188//
189
190
191template<class Scalar>
193 :dumpSensitivities_(false),
194 p_index_(-1),
195 g_index_(-1),
196 response_func_supports_x_dot_(false),
197 response_func_supports_D_x_dot_(false),
198 response_func_supports_D_p_(false)
199{}
200
201
202template<class Scalar>
204 const RCP<const Thyra::ModelEvaluator<Scalar> > &responseFunc,
205 const Thyra::ModelEvaluatorBase::InArgs<Scalar> &basePoint,
206 const int p_index,
207 const int g_index
208 )
209{
210
211 typedef Thyra::ModelEvaluatorBase MEB;
212
213 // ToDo: Validate input!
214
215 responseFunc_ = responseFunc;
216 basePoint_ = basePoint;
217 p_index_ = p_index;
218 g_index_ = g_index;
219
220 p_space_ = responseFunc_->get_p_space(p_index_);
221 g_space_ = responseFunc_->get_g_space(g_index_);
222
223
224 MEB::InArgs<Scalar>
225 responseInArgs = responseFunc_->createInArgs();
226 response_func_supports_x_dot_ =
227 responseInArgs.supports(MEB::IN_ARG_x_dot);
228 MEB::OutArgs<Scalar>
229 responseOutArgs = responseFunc_->createOutArgs();
230 response_func_supports_D_x_dot_ =
231 !responseOutArgs.supports(MEB::OUT_ARG_DgDx_dot,g_index_).none();
232 response_func_supports_D_p_ =
233 !responseOutArgs.supports(MEB::OUT_ARG_DgDp,g_index_,p_index_).none();
234
235 clearCache();
236
237}
238
239
240template<class Scalar>
242 const RCP<const Thyra::ModelEvaluator<Scalar> > &responseFunc,
243 const Thyra::ModelEvaluatorBase::InArgs<Scalar> &basePoint
244 )
245{
246 // ToDo: Validate that responseFunc is the same structure as the one already
247 // set!
248 responseFunc_ = responseFunc;
249 basePoint_ = basePoint;
250}
251
252
253template<class Scalar>
254const RCP<Thyra::VectorBase<Scalar> >
256{
257 return Thyra::createMember(g_space_);
258}
259
260
261template<class Scalar>
262const RCP<Thyra::MultiVectorBase<Scalar> >
264{
265 return Thyra::createMembers(g_space_,p_space_->dim());
266}
267
268
269template<class Scalar>
271 const Thyra::VectorBase<Scalar> *x_dot,
272 const Thyra::VectorBase<Scalar> &x,
273 const Scalar t,
274 Thyra::VectorBase<Scalar> *g_hat
275 ) const
276{
277 computeResponseAndSensitivityImpl(x_dot,0,x,0,t,g_hat,0);
278}
279
280
281template<class Scalar>
283 const Thyra::VectorBase<Scalar> *x_dot,
284 const Thyra::MultiVectorBase<Scalar> *S_dot,
285 const Thyra::VectorBase<Scalar> &x,
286 const Thyra::MultiVectorBase<Scalar> &S,
287 const Scalar t,
288 Thyra::VectorBase<Scalar> *g_hat,
289 Thyra::MultiVectorBase<Scalar> *D_g_hat_D_p
290 ) const
291{
292 computeResponseAndSensitivityImpl(x_dot,S_dot,x,&S,t,g_hat,D_g_hat_D_p);
293}
294
295
296// private
297
298
299template<class Scalar>
301{
302 D_g_D_x_dot_ = Teuchos::null;
303 D_g_D_x_ = Teuchos::null;
304 D_g_D_p_ = Teuchos::null;
305}
306
307
308template<class Scalar>
309void ForwardResponseSensitivityComputer<Scalar>::createCache(
310 const bool computeSens
311 ) const
312{
313 if (computeSens) {
314 if (response_func_supports_D_x_dot_ && is_null(D_g_D_x_dot_))
315 D_g_D_x_dot_ = responseFunc_->create_DgDx_dot_op(g_index_);
316 D_g_D_x_ = responseFunc_->create_DgDx_op(g_index_);
317 if (response_func_supports_D_p_ && is_null(D_g_D_p_))
318 D_g_D_p_ = createMembers(g_space_,p_space_->dim());
319 }
320}
321
322
323template<class Scalar>
324void ForwardResponseSensitivityComputer<Scalar>::computeResponseAndSensitivityImpl(
325 const Thyra::VectorBase<Scalar> *x_dot,
326 const Thyra::MultiVectorBase<Scalar> *S_dot,
327 const Thyra::VectorBase<Scalar> &x,
328 const Thyra::MultiVectorBase<Scalar> *S,
329 const Scalar t,
330 Thyra::VectorBase<Scalar> *g_hat,
331 Thyra::MultiVectorBase<Scalar> *D_g_hat_D_p
332 ) const
333{
334
335 using Teuchos::rcp;
336 using Teuchos::ptr;
337 using Teuchos::includesVerbLevel;
338 typedef ScalarTraits<Scalar> ST;
339 using Thyra::apply;
340 using Thyra::Vp_V;
341 typedef Thyra::ModelEvaluatorBase MEB;
342
343 //
344 // A) Setup for output
345 //
346
347 const RCP<Teuchos::FancyOStream> out = this->getOStream();
348 const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
349
350 const bool trace =
351 out.get() && includesVerbLevel(verbLevel,Teuchos::VERB_LOW);
352 const bool print_norms =
353 out.get() && includesVerbLevel(verbLevel,Teuchos::VERB_MEDIUM);
354 const bool print_x =
355 out.get() && includesVerbLevel(verbLevel,Teuchos::VERB_EXTREME);
356
357 //
358 // B) Initialize storage
359 //
360
361 const bool computeSens = ( D_g_hat_D_p != 0 );
362 createCache(computeSens);
363
364 //
365 // C) Evaluate the response function
366 //
367
368 //
369 // C.1) Setup input/output and evaluate the response function
370 //
371
372 if (trace)
373 *out << "\nEvaluating response function at time t = " << t << " ...\n";
374
375 // C.1.a) Setup the input arguments
376
377 MEB::InArgs<Scalar> responseInArgs = responseFunc_->createInArgs();
378 responseInArgs.setArgs(basePoint_);
379 responseInArgs.set_x(rcp(&x,false));
380 if (response_func_supports_x_dot_)
381 responseInArgs.set_x_dot(rcp(x_dot,false));
382
383 // C.1.b) Setup output arguments
384
385 MEB::OutArgs<Scalar> responseOutArgs = responseFunc_->createOutArgs();
386
387 if (g_hat)
388 responseOutArgs.set_g(g_index_,rcp(g_hat,false));
389
390 if (computeSens) {
391
392 // D_g_D_x_dot
393 if (response_func_supports_D_x_dot_) {
394 responseOutArgs.set_DgDx_dot(
395 g_index_,
396 MEB::Derivative<Scalar>(D_g_D_x_dot_)
397 );
398 }
399
400 // D_g_D_x
401 responseOutArgs.set_DgDx(
402 g_index_,
403 MEB::Derivative<Scalar>(D_g_D_x_)
404 );
405
406 // D_g_D_p
407 if (response_func_supports_D_p_) {
408 responseOutArgs.set_DgDp(
409 g_index_, p_index_,
410 MEB::Derivative<Scalar>(D_g_D_p_,MEB::DERIV_MV_BY_COL)
411 );
412 }
413
414 }
415
416 // C.1.c) Evaluate the response function k
417
418 {
419 RYTHMOS_FUNC_TIME_MONITOR("Rythmos:ForwardResponseSensitivityComputer::evalModel: evalResponse");
420 responseFunc_->evalModel( responseInArgs, responseOutArgs );
421 }
422
423 // C.1.d) Print the outputs just coputed
424
425 if (print_norms) {
426 if (g_hat)
427 *out << "\n||g_hat||inf = " << norm_inf(*g_hat) << std::endl;
428 if (computeSens && !is_null(D_g_D_p_))
429 *out << "\n||D_g_D_p_||inf = " << norms_inf(*D_g_D_p_) << std::endl;
430 }
431
432 if ( g_hat && (dumpSensitivities_ || print_x) )
433 *out << "\ng_hat = " << *g_hat;
434
435 if (computeSens && print_x) {
436 if (!is_null(D_g_D_x_dot_))
437 *out << "\nD_g_D_x_dot = " << *D_g_D_x_ << std::endl;
438 if (!is_null(D_g_D_x_))
439 *out << "\nD_g_D_x = " << *D_g_D_x_ << std::endl;
440 if (!is_null(D_g_D_p_))
441 *out << "\nD_g_D_p = " << *D_g_D_p_ << std::endl;
442 }
443
444 //
445 // C.2) Assemble the output response function sensitivity D_d_hat_D_p
446 //
447
448 // D_g_hat_D_p = DgDx_dot * S_dot + DgDx * S + DgDp
449
450 if (computeSens) {
451
452 RYTHMOS_FUNC_TIME_MONITOR("Rythmos:ForwardResponseSensitivityComputer::evalModel: computeSens");
453
454 if (trace)
455 *out << "\nD_g_hat_D_p = DgDx_dot * S_dot + DgDx * S + DgDp ...\n";
456
457 assign( ptr(D_g_hat_D_p), ST::zero() );
458
459 // D_g_hat_D_p += DgDx_dot * S_dot
460 if (response_func_supports_D_x_dot_) {
461 apply( *D_g_D_x_dot_, Thyra::NOTRANS, *S_dot,
462 ptr(D_g_hat_D_p), ST::one(), ST::one() );
463 }
464
465 // D_g_hat_D_p += DgDx * S
466 apply( *D_g_D_x_, Thyra::NOTRANS, *S,
467 ptr(D_g_hat_D_p), ST::one(), ST::one() );
468
469 // D_g_hat_D_p += DgDp
470 if (response_func_supports_D_p_) {
471 Vp_V( ptr(D_g_hat_D_p), *D_g_D_p_ );
472 }
473
474 if (dumpSensitivities_ || print_x)
475 *out << "\nD_g_hat_D_p = "
476 << Teuchos::describe(*D_g_hat_D_p,Teuchos::VERB_EXTREME);
477
478 }
479
480}
481
482
483} // namespace Rythmos
484
485
486#endif // RYTHMOS_FORWARD_RESPONSE_SENSITIVITY_COMPUTER_HPP
Concrete utility class for computing (assembling) forward transient response sensitivities.
const RCP< Thyra::MultiVectorBase< Scalar > > create_D_g_hat_D_p() const
void computeResponseAndSensitivity(const Thyra::VectorBase< Scalar > *x_dot, const Thyra::MultiVectorBase< Scalar > *S_dot, const Thyra::VectorBase< Scalar > &x, const Thyra::MultiVectorBase< Scalar > &S, const Scalar t, Thyra::VectorBase< Scalar > *g_hat, Thyra::MultiVectorBase< Scalar > *D_g_hat_D_p) const
Compute the reduced sensitivity and perhaps the response itself at a point (xdot,x,...
void setResponseFunction(const RCP< const Thyra::ModelEvaluator< Scalar > > &responseFunc, const Thyra::ModelEvaluatorBase::InArgs< Scalar > &basePoint, const int p_index, const int g_index)
Set the response function for the first time.
const RCP< Thyra::VectorBase< Scalar > > create_g_hat() const
STANDARD_MEMBER_COMPOSITION_MEMBERS(bool, dumpSensitivities)
void computeResponse(const Thyra::VectorBase< Scalar > *x_dot, const Thyra::VectorBase< Scalar > &x, const Scalar t, Thyra::VectorBase< Scalar > *g_hat) const
Compute the reduced response at a point (xdot,x,t).
void resetResponseFunction(const RCP< const Thyra::ModelEvaluator< Scalar > > &responseFunc, const Thyra::ModelEvaluatorBase::InArgs< Scalar > &basePoint)
Reset the point-specific response function along with its base point.