Rythmos - Transient Integration for Differential Equations Version of the Day
Loading...
Searching...
No Matches
Rythmos_ForwardSensitivityExplicitModelEvaluator.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_SENSITIVITY_EXPLICIT_MODEL_EVALUATOR_HPP
30#define RYTHMOS_FORWARD_SENSITIVITY_EXPLICIT_MODEL_EVALUATOR_HPP
31
32
33#include "Rythmos_IntegratorBase.hpp"
34#include "Rythmos_ForwardSensitivityModelEvaluatorBase.hpp"
35#include "Thyra_ModelEvaluator.hpp" // Interface
36#include "Thyra_StateFuncModelEvaluatorBase.hpp" // Implementation
37#include "Thyra_ModelEvaluatorDelegatorBase.hpp"
38#include "Thyra_DefaultMultiVectorProductVectorSpace.hpp"
39#include "Thyra_DefaultMultiVectorProductVector.hpp"
40#include "Thyra_DefaultMultiVectorLinearOpWithSolve.hpp"
41#include "Thyra_VectorStdOps.hpp"
42#include "Thyra_MultiVectorStdOps.hpp"
43
44
45namespace Rythmos {
46
47
125template<class Scalar>
127 : virtual public ForwardSensitivityModelEvaluatorBase<Scalar>
128{
129public:
130
133
136
138
141
144 const RCP<const Thyra::ModelEvaluator<Scalar> > &stateModel,
145 const int p_index
146 );
147
150 const RCP<const Thyra::ModelEvaluator<Scalar> >& stateModel,
151 const RCP<const Thyra::VectorSpaceBase<Scalar> >& p_space
152 );
153
155 RCP<const Thyra::ModelEvaluator<Scalar> >
156 getStateModel() const;
157
159 RCP<Thyra::ModelEvaluator<Scalar> >
160 getNonconstStateModel() const;
161
163 int get_p_index() const;
164
166 RCP<const Thyra::DefaultMultiVectorProductVectorSpace<Scalar> >
167 get_s_bar_space() const;
168
170 RCP<const Thyra::VectorSpaceBase<Scalar> > get_p_sens_space() const;
171
175 Ptr<StepperBase<Scalar> > stateStepper,
176 bool forceUpToDateW
177 );
178
180
183
185 RCP<const Thyra::VectorSpaceBase<Scalar> > get_x_space() const;
187 RCP<const Thyra::VectorSpaceBase<Scalar> > get_f_space() const;
189 Thyra::ModelEvaluatorBase::InArgs<Scalar> getNominalValues() const;
191 RCP<Thyra::LinearOpWithSolveBase<Scalar> > create_W() const;
193 Thyra::ModelEvaluatorBase::InArgs<Scalar> createInArgs() const;
194
196
197private:
198
201
203 Thyra::ModelEvaluatorBase::OutArgs<Scalar> createOutArgsImpl() const;
205 void evalModelImpl(
206 const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs,
207 const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs
208 ) const;
209
211
212private:
213
214 // /////////////////////////
215 // Private data members
216
217 RCP<const Thyra::ModelEvaluator<Scalar> > stateModel_;
218 int p_index_;
219 int np_;
220
221 RCP<const Thyra::DefaultMultiVectorProductVectorSpace<Scalar> > s_bar_space_;
222 RCP<const Thyra::DefaultMultiVectorProductVectorSpace<Scalar> > f_sens_space_;
223
224 Thyra::ModelEvaluatorBase::InArgs<Scalar> nominalValues_;
225
226 mutable Thyra::ModelEvaluatorBase::InArgs<Scalar> stateBasePoint_;
227
228 mutable RCP<Thyra::LinearOpBase<Scalar> > DfDx_;
229 mutable RCP<Thyra::MultiVectorBase<Scalar> > DfDp_;
230
231 mutable RCP<Thyra::LinearOpBase<Scalar> > DfDx_compute_;
232 mutable RCP<Thyra::MultiVectorBase<Scalar> > DfDp_compute_;
233
234 // /////////////////////////
235 // Private member functions
236
237 void wrapNominalValuesAndBounds();
238
239 void computeDerivativeMatrices(
240 const Thyra::ModelEvaluatorBase::InArgs<Scalar> &point
241 ) const;
242
243};
244
245
246// /////////////////////////////////
247// Implementations
248
249
250// Constructors/Intializers/Accessors
251
252
253template<class Scalar>
254RCP<ForwardSensitivityExplicitModelEvaluator<Scalar> >
255forwardSensitivityExplicitModelEvaluator()
256{
257 RCP<ForwardSensitivityExplicitModelEvaluator<Scalar> > fseme =
259 return fseme;
260}
261
262
263template<class Scalar>
265 : p_index_(0), np_(-1)
266{}
267
268
269
270// Public functions overridden from ForwardSensitivityModelEvaluatorBase
271
272
273template<class Scalar>
275 const RCP<const Thyra::ModelEvaluator<Scalar> > &stateModel,
276 const int p_index
277 )
278{
279
280 typedef Thyra::ModelEvaluatorBase MEB;
281
282 //
283 // Validate input
284 //
285
286 TEUCHOS_TEST_FOR_EXCEPT( is_null(stateModel) );
287 TEUCHOS_TEST_FOR_EXCEPTION(
288 !( 0 <= p_index && p_index < stateModel->Np() ), std::logic_error,
289 "Error, p_index does not fall in the range [0,"<<(stateModel->Np()-1)<<"]!" );
290 // ToDo: Validate support for DfDp!
291
292 //
293 // Set the input objects
294 //
295
296 stateModel_ = stateModel;
297 p_index_ = p_index;
298 np_ = stateModel_->get_p_space(p_index)->dim();
299
300 //
301 // Create the structure of the model
302 //
303
304 s_bar_space_ = Thyra::multiVectorProductVectorSpace(
305 stateModel_->get_x_space(), np_
306 );
307
308 f_sens_space_ = Thyra::multiVectorProductVectorSpace(
309 stateModel_->get_f_space(), np_
310 );
311
312 nominalValues_ = this->createInArgs();
313
314 this->wrapNominalValuesAndBounds();
315
316 //
317 // Wipe out matrix storage
318 //
319
320 stateBasePoint_ = MEB::InArgs<Scalar>();
321 DfDx_ = Teuchos::null;
322 DfDp_ = Teuchos::null;
323 DfDx_compute_ = Teuchos::null;
324 DfDp_compute_ = Teuchos::null;
325
326}
327
328
329template<class Scalar>
331 const RCP<const Thyra::ModelEvaluator<Scalar> >& /* stateModel */,
332 const RCP<const Thyra::VectorSpaceBase<Scalar> >& /* p_space */
333 )
334{
335 TEUCHOS_TEST_FOR_EXCEPT_MSG(true, "ToDo: Implement initializeStructureInitCondOnly()!" );
336}
337
338
339template<class Scalar>
340RCP<const Thyra::ModelEvaluator<Scalar> >
342{
343 return stateModel_;
344}
345
346
347template<class Scalar>
348RCP<Thyra::ModelEvaluator<Scalar> >
350{
351 return Teuchos::null;
352}
353
354
355template<class Scalar>
357{
358 return p_index_;
359}
360
361
362template<class Scalar>
363RCP<const Thyra::DefaultMultiVectorProductVectorSpace<Scalar> >
365{
366 return s_bar_space_;
367}
368
369
370template<class Scalar>
371RCP<const Thyra::VectorSpaceBase<Scalar> >
373{
374 return stateModel_->get_p_space(p_index_);
375}
376
377
378template<class Scalar>
380 Ptr<StepperBase<Scalar> > stateStepper,
381 bool /* forceUpToDateW */
382 )
383{
384 TEUCHOS_ASSERT( Teuchos::nonnull(stateStepper) );
385#ifdef HAVE_RYTHMOS_DEBUG
386 TEUCHOS_TEST_FOR_EXCEPTION(
387 is_null(stateModel_), std::logic_error,
388 "Error, you must call intializeStructure(...) before you call initializePointState(...)"
389 );
390#endif // HAVE_RYTHMOS_DEBUG
391
392 Scalar curr_t = stateStepper->getStepStatus().time;
393 RCP<const Thyra::VectorBase<Scalar> > x;
394 x = get_x(*stateStepper,curr_t);
395#ifdef HAVE_RYTHMOS_DEBUG
396 TEUCHOS_TEST_FOR_EXCEPT( Teuchos::is_null(x) );
397#endif // HAVE_RYTHMOS_DEBUG
398
399 stateBasePoint_ = stateStepper->getInitialCondition(); // set parameters
400 stateBasePoint_.set_x( x );
401 stateBasePoint_.set_t( curr_t );
402
403 // Set whatever derivatives where passed in. If an input in null, then the
404 // member will be null and the null linear operators will be computed later
405 // just in time.
406
407 wrapNominalValuesAndBounds();
408
409}
410
411
412// Public functions overridden from ModelEvaulator
413
414
415template<class Scalar>
416RCP<const Thyra::VectorSpaceBase<Scalar> >
418{
419 return s_bar_space_;
420}
421
422
423template<class Scalar>
424RCP<const Thyra::VectorSpaceBase<Scalar> >
426{
427 return f_sens_space_;
428}
429
430
431template<class Scalar>
432Thyra::ModelEvaluatorBase::InArgs<Scalar>
434{
435 return nominalValues_;
436}
437
438
439template<class Scalar>
440RCP<Thyra::LinearOpWithSolveBase<Scalar> >
442{
443 return Thyra::multiVectorLinearOpWithSolve<Scalar>();
444}
445
446
447template<class Scalar>
448Thyra::ModelEvaluatorBase::InArgs<Scalar>
450{
451 TEUCHOS_ASSERT( !is_null(stateModel_) );
452 typedef Thyra::ModelEvaluatorBase MEB;
453 MEB::InArgs<Scalar> stateModelInArgs = stateModel_->createInArgs();
454 MEB::InArgsSetup<Scalar> inArgs;
455 inArgs.setModelEvalDescription(this->description());
456 inArgs.setSupports( MEB::IN_ARG_x );
457 inArgs.setSupports( MEB::IN_ARG_t );
458 inArgs.setSupports( MEB::IN_ARG_beta,
459 stateModelInArgs.supports(MEB::IN_ARG_beta) );
460 return inArgs;
461}
462
463
464// Private functions overridden from ModelEvaulatorDefaultBase
465
466
467template<class Scalar>
468Thyra::ModelEvaluatorBase::OutArgs<Scalar>
470{
471 TEUCHOS_ASSERT( !is_null(stateModel_) );
472 typedef Thyra::ModelEvaluatorBase MEB;
473
474 MEB::OutArgs<Scalar> stateModelOutArgs = stateModel_->createOutArgs();
475 MEB::OutArgsSetup<Scalar> outArgs;
476
477 outArgs.setModelEvalDescription(this->description());
478
479 outArgs.setSupports(MEB::OUT_ARG_f);
480
481 return outArgs;
482
483}
484
485
486template<class Scalar>
487void ForwardSensitivityExplicitModelEvaluator<Scalar>::evalModelImpl(
488 const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs,
489 const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs
490 ) const
491{
492
493 using Teuchos::rcp_dynamic_cast;
494 typedef Teuchos::ScalarTraits<Scalar> ST;
495 // typedef Thyra::ModelEvaluatorBase MEB; // unused
496 typedef Teuchos::VerboseObjectTempState<Thyra::ModelEvaluatorBase> VOTSME;
497
498 THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_GEN_BEGIN(
499 "ForwardSensitivityExplicitModelEvaluator", inArgs, outArgs, Teuchos::null );
500
501 //
502 // Update the derivative matrices if they are not already updated for the
503 // given time!.
504 //
505
506 {
507 RYTHMOS_FUNC_TIME_MONITOR_DIFF(
508 "Rythmos:ForwardSensitivityExplicitModelEvaluator::evalModel: computeMatrices",
509 Rythmos_FSEME
510 );
511 computeDerivativeMatrices(inArgs);
512 }
513
514 //
515 // InArgs
516 //
517
518 RCP<const Thyra::DefaultMultiVectorProductVector<Scalar> >
519 s_bar = rcp_dynamic_cast<const Thyra::DefaultMultiVectorProductVector<Scalar> >(
520 inArgs.get_x().assert_not_null(), true
521 );
522 RCP<const Thyra::MultiVectorBase<Scalar> >
523 S = s_bar->getMultiVector();
524
525 //
526 // OutArgs
527 //
528
529 RCP<Thyra::DefaultMultiVectorProductVector<Scalar> >
530 f_sens = rcp_dynamic_cast<Thyra::DefaultMultiVectorProductVector<Scalar> >(
531 outArgs.get_f(), true
532 );
533
534 RCP<Thyra::MultiVectorBase<Scalar> >
535 F_sens = f_sens->getNonconstMultiVector().assert_not_null();
536
537 //
538 // Compute the requested functions
539 //
540
541 if(!is_null(F_sens)) {
542
543 RYTHMOS_FUNC_TIME_MONITOR_DIFF(
544 "Rythmos:ForwardSensitivityExplicitModelEvaluator::evalModel: computeSens",
545 Rythmos_FSEME
546 );
547
548 // Form the residual: df/dx * S + df/dp
549 // F_sens = df/dx * S
550 Thyra::apply(
551 *DfDx_, Thyra::NOTRANS,
552 *S, F_sens.ptr(),
553 ST::one(), ST::zero()
554 );
555 // F_sens += d(f)/d(p)
556 Vp_V( F_sens.ptr(), *DfDp_ );
557 }
558
559 THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_END();
560
561}
562
563
564// private
565
566
567template<class Scalar>
568void ForwardSensitivityExplicitModelEvaluator<Scalar>::wrapNominalValuesAndBounds()
569{
570 TEUCHOS_ASSERT( !is_null(stateModel_) );
571
572 // using Teuchos::rcp_dynamic_cast; // unused
573 // typedef Thyra::ModelEvaluatorBase MEB; // unused
574 typedef Teuchos::ScalarTraits<Scalar> ST;
575
576 // nominalValues_.clear(); // ToDo: Implement this!
577
578 nominalValues_.set_t(stateModel_->getNominalValues().get_t());
579
580 // 2007/05/22: rabartl: Note: Currently there is not much of a reason to set
581 // an initial condition here since the initial condition for the
582 // sensitivities is really being set in the ForwardSensitivityStepper
583 // object! In the future, a more general use of this class might benefit
584 // from setting the initial condition here.
585
586 // 2009/07/20: tscoffe/ccober: This is the future. We're going to use this
587 // in a more general way, so we need legitimate nominal values now.
588 RCP<VectorBase<Scalar> > s_bar_ic = Thyra::createMember(this->get_x_space());
589 Thyra::V_S(s_bar_ic.ptr(),ST::zero());
590 nominalValues_.set_x(s_bar_ic);
591}
592
593
594template<class Scalar>
595void ForwardSensitivityExplicitModelEvaluator<Scalar>::computeDerivativeMatrices(
596 const Thyra::ModelEvaluatorBase::InArgs<Scalar> &/* point */
597 ) const
598{
599 TEUCHOS_ASSERT( !is_null(stateModel_) );
600
601 typedef Thyra::ModelEvaluatorBase MEB;
602 typedef Teuchos::VerboseObjectTempState<MEB> VOTSME;
603
604 Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
605 Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
606
607 MEB::InArgs<Scalar> inArgs = stateBasePoint_;
608 MEB::OutArgs<Scalar> outArgs = stateModel_->createOutArgs();
609
610 if (is_null(DfDx_)) {
611 DfDx_ = stateModel_->create_W_op();
612 }
613 if (inArgs.supports(MEB::IN_ARG_beta)) {
614 inArgs.set_beta(1.0);
615 }
616 outArgs.set_W_op(DfDx_);
617
618 if (is_null(DfDp_)) {
619 DfDp_ = Thyra::create_DfDp_mv(
620 *stateModel_,p_index_,
621 MEB::DERIV_MV_BY_COL
622 ).getMultiVector();
623 }
624 outArgs.set_DfDp(
625 p_index_,
626 MEB::Derivative<Scalar>(DfDp_,MEB::DERIV_MV_BY_COL)
627 );
628
629 VOTSME stateModel_outputTempState(stateModel_,out,verbLevel);
630 stateModel_->evalModel(inArgs,outArgs);
631
632
633}
634
635
636} // namespace Rythmos
637
638
639#endif // RYTHMOS_FORWARD_SENSITIVITY_EXPLICIT_MODEL_EVALUATOR_HPP
Explicit forward sensitivity transient ModelEvaluator subclass.
RCP< const Thyra::VectorSpaceBase< Scalar > > get_p_sens_space() const
void initializeStructureInitCondOnly(const RCP< const Thyra::ModelEvaluator< Scalar > > &stateModel, const RCP< const Thyra::VectorSpaceBase< Scalar > > &p_space)
void initializeStructure(const RCP< const Thyra::ModelEvaluator< Scalar > > &stateModel, const int p_index)
RCP< const Thyra::DefaultMultiVectorProductVectorSpace< Scalar > > get_s_bar_space() const
void initializePointState(Ptr< StepperBase< Scalar > > stateStepper, bool forceUpToDateW)
Initialize full state for a single point in time.
Forward sensitivity transient ModelEvaluator node interface class.
Base class for defining stepper functionality.