Rythmos - Transient Integration for Differential Equations Version of the Day
Loading...
Searching...
No Matches
Rythmos_ForwardSensitivityIntegratorAsModelEvaluator.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
30#ifndef RYTHMOS_FORWARD_SENSITIVITY_INTEGRATOR_AS_MODEL_EVALUATOR_HPP
31#define RYTHMOS_FORWARD_SENSITIVITY_INTEGRATOR_AS_MODEL_EVALUATOR_HPP
32
33
34#include "Thyra_ResponseOnlyModelEvaluatorBase.hpp"
35
36#include "Rythmos_StepperBase.hpp"
37#include "Rythmos_IntegratorBase.hpp"
38#include "Rythmos_InterpolationBufferHelpers.hpp"
39#include "Rythmos_ForwardSensitivityStepper.hpp"
40#include "Rythmos_ForwardResponseSensitivityComputer.hpp"
41#include "Rythmos_extractStateAndSens.hpp"
42#include "Thyra_ModelEvaluatorDelegatorBase.hpp"
43#include "Thyra_DefaultMultiVectorProductVectorSpace.hpp"
44#include "Thyra_DefaultProductVectorSpace.hpp"
45#include "Teuchos_ParameterListAcceptor.hpp"
46#include "Teuchos_VerboseObjectParameterListHelpers.hpp"
47#include "Teuchos_Assert.hpp"
48
49
50namespace Rythmos {
51
52
53namespace ForwardSensitivityIntegratorAsModelEvaluatorTypes {
55enum EResponseType {
57 RESPONSE_TYPE_SUM,
59 RESPONSE_TYPE_BLOCK
60};
61} // namespace ForwardSensitivityIntegratorAsModelEvaluatorTypes
62
63
138template<class Scalar>
140 : virtual public Thyra::ResponseOnlyModelEvaluatorBase<Scalar>
141 , virtual public Teuchos::ParameterListAcceptor
142{
143public:
144
146
149
152
206 void initialize(
207 const RCP<StepperBase<Scalar> > &stateStepper,
208 const RCP<IntegratorBase<Scalar> > &stateIntegrator,
209 const RCP<ForwardSensitivityStepper<Scalar> > &stateAndSensStepper,
210 const RCP<IntegratorBase<Scalar> > &stateAndSensIntegrator,
211 const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateAndSensInitCond,
212 const Array<Scalar> &responseTimes,
213 const Array<RCP<const Thyra::ModelEvaluator<Scalar> > > &responseFuncs,
214 const Array<Thyra::ModelEvaluatorBase::InArgs<Scalar> > &responseFuncBasePoints,
215 const ForwardSensitivityIntegratorAsModelEvaluatorTypes::EResponseType responseType
216 );
217
219 const Array<Scalar>& getResponseTimes() const;
220
222
223
226
234 void setParameterList(RCP<Teuchos::ParameterList> const& paramList);
236 RCP<Teuchos::ParameterList> getNonconstParameterList();
238 RCP<Teuchos::ParameterList> unsetParameterList();
240 RCP<const Teuchos::ParameterList> getParameterList() const;
249 RCP<const Teuchos::ParameterList> getValidParameters() const;
250
252
255
257 int Np() const;
259 int Ng() const;
261 RCP<const Thyra::VectorSpaceBase<Scalar> > get_p_space(int l) const;
263 RCP<const Thyra::VectorSpaceBase<Scalar> > get_g_space(int j) const;
265 Thyra::ModelEvaluatorBase::InArgs<Scalar> createInArgs() const;
266
268
269private:
270
273
275 Thyra::ModelEvaluatorBase::OutArgs<Scalar> createOutArgsImpl() const;
277 void evalModelImpl(
278 const Thyra::ModelEvaluatorBase::InArgs<Scalar>& inArgs,
279 const Thyra::ModelEvaluatorBase::OutArgs<Scalar>& outArgs
280 ) const;
281
283
284private:
285
286 // //////////////////////
287 // Private types
288
289 typedef Teuchos::Array<RCP<const Thyra::VectorSpaceBase<Scalar> > > SpaceArray_t;
290
291
292 // //////////////////////
293 // Private data members
294
295 RCP<Teuchos::ParameterList> paramList_;
296 bool dumpSensitivities_;
297
298 RCP<StepperBase<Scalar> > stateStepper_;
299 RCP<IntegratorBase<Scalar> > stateIntegrator_;
300 mutable Thyra::ModelEvaluatorBase::InArgs<Scalar> stateInitCond_;
301
302 RCP<ForwardSensitivityStepper<Scalar> > stateAndSensStepper_;
303 RCP<IntegratorBase<Scalar> > stateAndSensIntegrator_;
304 mutable Thyra::ModelEvaluatorBase::InArgs<Scalar> stateAndSensInitCond_;
305
306 Array<Scalar> responseTimes_;
307 Array<RCP<const Thyra::ModelEvaluator<Scalar> > > responseFuncs_;
308 mutable Array<Thyra::ModelEvaluatorBase::InArgs<Scalar> > responseFuncInArgs_;
309 ForwardSensitivityIntegratorAsModelEvaluatorTypes::EResponseType responseType_;
310 bool response_func_supports_x_dot_;
311 bool response_func_supports_D_x_dot_;
312 bool response_func_supports_D_p_;
313
314 int p_index_;
315 int g_index_;
316 Scalar finalTime_;
317
318 int Np_;
319 int Ng_;
320
321 SpaceArray_t g_space_;
322 SpaceArray_t p_space_;
323
324 static const std::string dumpSensitivities_name_;
325 static const bool dumpSensitivities_default_;
326
327};
328
329
331template<class Scalar>
332RCP<ForwardSensitivityIntegratorAsModelEvaluator<Scalar> >
333forwardSensitivityIntegratorAsModelEvaluator(
334 const RCP<StepperBase<Scalar> > &stateStepper,
335 const RCP<IntegratorBase<Scalar> > &stateIntegrator,
336 const RCP<ForwardSensitivityStepper<Scalar> > &stateAndSensStepper,
337 const RCP<IntegratorBase<Scalar> > &stateAndSensIntegrator,
338 const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateAndSensInitCond,
339 const Array<Scalar> &responseTimes,
340 const Array<RCP<const Thyra::ModelEvaluator<Scalar> > > &responseFuncs,
341 const Array<Thyra::ModelEvaluatorBase::InArgs<Scalar> > &responseFuncBasePoints,
342 const ForwardSensitivityIntegratorAsModelEvaluatorTypes::EResponseType responseType
343 )
344{
345 using Teuchos::rcp;
346 RCP<ForwardSensitivityIntegratorAsModelEvaluator<Scalar> >
347 sensIntegratorAsModelEvaluator = rcp(new ForwardSensitivityIntegratorAsModelEvaluator<Scalar>());
348 sensIntegratorAsModelEvaluator->initialize(
349 stateStepper, stateIntegrator,
350 stateAndSensStepper, stateAndSensIntegrator,
351 stateAndSensInitCond,
352 responseTimes, responseFuncs,
353 responseFuncBasePoints, responseType
354 );
355 return sensIntegratorAsModelEvaluator;
356}
357
358
359// ////////////////////////
360// Definitions
361
362
363// Static members
364
365
366template<class Scalar>
367const std::string
368ForwardSensitivityIntegratorAsModelEvaluator<Scalar>::dumpSensitivities_name_
369= "Dump Sensitivities";
370
371template<class Scalar>
372const bool
373ForwardSensitivityIntegratorAsModelEvaluator<Scalar>::dumpSensitivities_default_
374= false;
375
376
377// Constructors, Initialization, Misc.
378
379
380template<class Scalar>
382 :dumpSensitivities_(dumpSensitivities_default_),
383 responseType_(ForwardSensitivityIntegratorAsModelEvaluatorTypes::RESPONSE_TYPE_SUM),
384 response_func_supports_x_dot_(false),
385 response_func_supports_D_x_dot_(false),
386 response_func_supports_D_p_(false),
387 p_index_(-1),
388 g_index_(-1),
389 finalTime_(-std::numeric_limits<Scalar>::max()),
390 Np_(0),
391 Ng_(0)
392{}
393
394
395template<class Scalar>
397 const RCP<StepperBase<Scalar> > &stateStepper,
398 const RCP<IntegratorBase<Scalar> > &stateIntegrator,
399 const RCP<ForwardSensitivityStepper<Scalar> > &stateAndSensStepper,
400 const RCP<IntegratorBase<Scalar> > &stateAndSensIntegrator,
401 const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateAndSensInitCond,
402 const Array<Scalar> &responseTimes,
403 const Array<RCP<const Thyra::ModelEvaluator<Scalar> > > &responseFuncs,
404 const Array<Thyra::ModelEvaluatorBase::InArgs<Scalar> > &responseFuncBasePoints,
405 const ForwardSensitivityIntegratorAsModelEvaluatorTypes::EResponseType responseType
406 )
407{
408
409 using Teuchos::as;
410 typedef Thyra::ModelEvaluatorBase MEB;
411 namespace FSIAMET = ForwardSensitivityIntegratorAsModelEvaluatorTypes;
412
413 //
414 // A) Validate and set input
415 //
416
417#ifdef HAVE_RYTHMOS_DEBUG
418 const int numResponseTimes = responseTimes.size();
419
420 TEUCHOS_TEST_FOR_EXCEPT(is_null(stateStepper));
421 TEUCHOS_TEST_FOR_EXCEPT(is_null(stateIntegrator));
422 TEUCHOS_TEST_FOR_EXCEPT(is_null(stateAndSensStepper));
423 TEUCHOS_TEST_FOR_EXCEPT(is_null(stateAndSensInitCond.get_x()));
424 TEUCHOS_TEST_FOR_EXCEPT(is_null(stateAndSensInitCond.get_x_dot()));
425 TEUCHOS_TEST_FOR_EXCEPT( !( numResponseTimes > 0 ) );
426 assertTimePointsAreSorted(responseTimes);
427 TEUCHOS_TEST_FOR_EXCEPT( as<int>(responseFuncs.size()) != numResponseTimes );
428 TEUCHOS_TEST_FOR_EXCEPT( as<int>(responseFuncBasePoints.size()) != numResponseTimes );
429 // ToDo: Assert that all of the observation models have the same response
430 // function spaces so that they can be added together!
431#endif // HAVE_RYTHMOS_DEBUG
432
433 stateStepper_ = stateStepper;
434 stateAndSensStepper_ = stateAndSensStepper;
435 stateAndSensInitCond_ = stateAndSensInitCond;
436 stateIntegrator_ = stateIntegrator;
437 if (!is_null(stateAndSensIntegrator))
438 stateAndSensIntegrator_ = stateAndSensIntegrator;
439 else
440 stateAndSensIntegrator_ = stateIntegrator_->cloneIntegrator().assert_not_null();
441 responseTimes_ = responseTimes;
442 responseFuncs_ = responseFuncs;
443 responseFuncInArgs_ = responseFuncBasePoints;
444 responseType_ = responseType;
445
446 finalTime_ = responseTimes_.back();
447
448 //
449 // B) Set the initial condition for the state-only problem
450 //
451
452 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >
453 stateModel = stateStepper_->getModel();
454
455 // Get base point pased in with no x or x_dot
456 MEB::InArgs<Scalar>
457 basePoint_no_x = stateAndSensInitCond_;
458 basePoint_no_x.set_x(Teuchos::null);
459 basePoint_no_x.set_x_dot(Teuchos::null);
460
461 // Create an empty InArgs for the state model/stepper
462 stateInitCond_ = stateModel->createInArgs();
463
464 // Set the base point (except x, x_dot).
465 stateInitCond_.setArgs(basePoint_no_x);
466
467 // Set x and x_dot
468 const RCP<const Thyra::ProductVectorBase<Scalar> >
469 x_bar_init = Thyra::productVectorBase<Scalar>(
470 stateAndSensInitCond_.get_x()
471 ),
472 x_bar_dot_init = Thyra::productVectorBase<Scalar>(
473 stateAndSensInitCond_.get_x_dot()
474 );
475 stateInitCond_.set_x(x_bar_init->getVectorBlock(0));
476 stateInitCond_.set_x_dot(x_bar_dot_init->getVectorBlock(0));
477
478 //
479 // C) Set up the info for this model evaluators interface
480 //
481
482 Np_ = 1;
483 p_index_ = getParameterIndex(*stateAndSensStepper_);
484 p_space_.clear();
485 p_space_.push_back(stateModel->get_p_space(p_index_));
486
487 Ng_ = 1;
488 g_index_ = 0; // ToDo: Accept this from input!
489 g_space_.clear();
490
491 if (responseType_ == FSIAMET::RESPONSE_TYPE_SUM) {
492 g_space_.push_back(responseFuncs[0]->get_g_space(0));
493 }
494 else if (responseType_ == FSIAMET::RESPONSE_TYPE_BLOCK) {
495 g_space_.push_back(
496 Thyra::productVectorSpace(
497 responseFuncs[0]->get_g_space(g_index_), responseFuncs.size()
498 )
499 );
500 }
501
502 MEB::InArgs<Scalar>
503 responseInArgs = responseFuncs[0]->createInArgs();
504 response_func_supports_x_dot_ =
505 responseInArgs.supports(MEB::IN_ARG_x_dot);
506 MEB::OutArgs<Scalar>
507 responseOutArgs = responseFuncs[0]->createOutArgs();
508 response_func_supports_D_x_dot_ =
509 !responseOutArgs.supports(MEB::OUT_ARG_DgDx_dot,g_index_).none();
510 response_func_supports_D_p_ =
511 !responseOutArgs.supports(MEB::OUT_ARG_DgDp,g_index_,p_index_).none();
512
513}
514
515
516template<class Scalar>
517const Array<Scalar>&
522
523
524// Overridden from Teuchos::ParameterListAcceptor
525
526
527template<class Scalar>
529 RCP<Teuchos::ParameterList> const& paramList
530 )
531{
532 TEUCHOS_TEST_FOR_EXCEPT(0==paramList.get());
533 paramList->validateParameters(*getValidParameters());
534 paramList_ = paramList;
535 dumpSensitivities_ = paramList_->get(
536 dumpSensitivities_name_, dumpSensitivities_default_);
537 Teuchos::readVerboseObjectSublist(&*paramList_,this);
538#ifdef HAVE_RYTHMOS_DEBUG
539 paramList_->validateParameters(*getValidParameters());
540#endif // HAVE_RYTHMOS_DEBUG
541
542}
543
544
545template<class Scalar>
546RCP<Teuchos::ParameterList>
551
552
553template<class Scalar>
554RCP<Teuchos::ParameterList>
556{
557 RCP<Teuchos::ParameterList> _paramList = paramList_;
558 paramList_ = Teuchos::null;
559 return _paramList;
560}
561
562
563template<class Scalar>
564RCP<const Teuchos::ParameterList>
569
570
571template<class Scalar>
572RCP<const Teuchos::ParameterList>
574{
575 static RCP<ParameterList> validPL;
576 if (is_null(validPL)) {
577 RCP<ParameterList> pl = Teuchos::parameterList();
578 pl->set(dumpSensitivities_name_, dumpSensitivities_default_,
579 "Set to true to have the individual sensitivities printed to\n"
580 "the output stream."
581 );
582 Teuchos::setupVerboseObjectSublist(&*pl);
583 validPL = pl;
584 }
585 return validPL;
586}
587
588
589// Public functions overridden from ModelEvaulator
590
591
592template<class Scalar>
597
598
599template<class Scalar>
604
605
606template<class Scalar>
607RCP<const Thyra::VectorSpaceBase<Scalar> >
609{
610#ifdef HAVE_RYTHMOS_DEBUG
611 TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE( l, 0, Np_ );
612#endif
613 return p_space_[l];
614}
615
616
617template<class Scalar>
618RCP<const Thyra::VectorSpaceBase<Scalar> >
620{
621#ifdef HAVE_RYTHMOS_DEBUG
622 TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE( j, 0, Ng_ );
623#endif
624 return g_space_[j];
625}
626
627
628template<class Scalar>
629Thyra::ModelEvaluatorBase::InArgs<Scalar>
631{
632 typedef Thyra::ModelEvaluatorBase MEB;
633 MEB::InArgsSetup<Scalar> inArgs;
634 inArgs.setModelEvalDescription(this->description());
635 inArgs.set_Np(Np_);
636 return inArgs;
637}
638
639
640// Private functions overridden from ModelEvaulatorDefaultBase
641
642
643template<class Scalar>
644Thyra::ModelEvaluatorBase::OutArgs<Scalar>
646{
647 typedef Thyra::ModelEvaluatorBase MEB;
648 typedef MEB::DerivativeSupport DS;
649 namespace FSIAMET = ForwardSensitivityIntegratorAsModelEvaluatorTypes;
650 MEB::OutArgsSetup<Scalar> outArgs;
651 outArgs.setModelEvalDescription(this->description());
652 outArgs.set_Np_Ng(Np_,Ng_);
653 outArgs.setSupports(MEB::OUT_ARG_DgDp, 0 ,0, MEB::DERIV_MV_BY_COL);
654 return outArgs;
655}
656
657
658template<class Scalar>
659void ForwardSensitivityIntegratorAsModelEvaluator<Scalar>::evalModelImpl(
660 const Thyra::ModelEvaluatorBase::InArgs<Scalar>& inArgs,
661 const Thyra::ModelEvaluatorBase::OutArgs<Scalar>& outArgs
662 ) const
663{
664
665 using Teuchos::as;
666 using Teuchos::null;
667 using Teuchos::describe;
668 using Teuchos::OSTab;
669 using Teuchos::rcp_dynamic_cast;
670 typedef Teuchos::ScalarTraits<Scalar> ST;
671 typedef Teuchos::VerboseObjectTempState<InterpolationBufferBase<Scalar> > VOTSSB;
672 typedef Thyra::ModelEvaluatorBase MEB;
673 namespace FSIAMET = ForwardSensitivityIntegratorAsModelEvaluatorTypes;
674
675 //
676 // Stream output and other basic setup
677 //
678
679 THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_GEN_BEGIN(
680 "Rythmos::ForwardSensitivityIntegratorAsModelEvaluator", inArgs, outArgs, null
681 );
682 VOTSSB stateIntegrator_outputTempState(
683 stateIntegrator_,out,incrVerbLevel(verbLevel,-1));
684 VOTSSB stateAndSensIntegrator_outputTempState(
685 stateAndSensIntegrator_,out,incrVerbLevel(verbLevel,-1));
686
687 const bool trace =
688 out.get() && includesVerbLevel(localVerbLevel,Teuchos::VERB_LOW);
689 const bool print_norms =
690 out.get() && includesVerbLevel(localVerbLevel,Teuchos::VERB_MEDIUM);
691 const bool print_x =
692 out.get() && includesVerbLevel(localVerbLevel,Teuchos::VERB_EXTREME);
693
694 const RCP<const Thyra::ModelEvaluator<Scalar> >
695 stateModel = stateStepper_->getModel();
696
697 //
698 // A) Process OutArgs first to see what functions we will be computing
699 //
700
701 RCP<Thyra::VectorBase<Scalar> >
702 d_hat = outArgs.get_g(0);
703
704 MEB::Derivative<Scalar>
705 D_d_hat_D_p = outArgs.get_DgDp(0,p_index_);
706
707 // Integrate state+sens or just state?
708 const bool integrateStateAndSens = !D_d_hat_D_p.isEmpty();
709
710 // Shortcut exit if no output was requested
711 if ( is_null(d_hat) && D_d_hat_D_p.isEmpty() ) {
712 if (trace)
713 *out << "\nSkipping evaluation since no outputs were requested!\n";
714 return;
715 }
716
717 //
718 // B) Process the InArgs knowing if we will integrate just the state or the
719 // state+sens.
720 //
721
722 const RCP<const Thyra::VectorBase<Scalar> >
723 p = inArgs.get_p(0).assert_not_null();
724
725 if (integrateStateAndSens) {
726 if (trace)
727 *out << "\nComputing D_d_hat_d_p by integrating state+sens ...\n";
728 stateAndSensInitCond_.set_p(p_index_,p);
729 }
730 else {
731 if (trace)
732 *out << "\nComputing just d_hat by integrating the state ...\n";
733 stateInitCond_.set_p(p_index_,p);
734 }
735
736 //
737 // C) Setup the stepper and the integrator for state+sens or only state
738 //
739
740 RCP<IntegratorBase<Scalar> > integrator;
741 if (integrateStateAndSens) {
742 stateAndSensStepper_->setInitialCondition(stateAndSensInitCond_);
743 stateAndSensIntegrator_->setStepper(stateAndSensStepper_,finalTime_);
744 integrator = stateAndSensIntegrator_;
745 }
746 else {
747 stateStepper_->setInitialCondition(stateInitCond_);
748 stateIntegrator_->setStepper(stateStepper_,finalTime_);
749 integrator = stateIntegrator_;
750 }
751
752 //
753 // D) Setup for computing and processing the individual response functions
754 //
755
756 ForwardResponseSensitivityComputer<Scalar>
757 forwardResponseSensitivityComputer;
758 forwardResponseSensitivityComputer.setOStream(out);
759 forwardResponseSensitivityComputer.setVerbLevel(localVerbLevel);
760 forwardResponseSensitivityComputer.dumpSensitivities(dumpSensitivities_);
761 forwardResponseSensitivityComputer.setResponseFunction(
762 responseFuncs_[0],
763 responseFuncs_[0]->createInArgs(), // Will replace this for each point!
764 p_index_, g_index_
765 );
766
767 // D.a) Create storage for the individual response function ouptuts g_k
768 // and its derivaitves
769
770 RCP<Thyra::VectorBase<Scalar> > g_k;
771 RCP<Thyra::MultiVectorBase<Scalar> > D_g_hat_D_p_k;
772
773 if (!is_null(d_hat)) {
774 g_k = forwardResponseSensitivityComputer.create_g_hat();
775 }
776 if (!D_d_hat_D_p.isEmpty()) {
777 D_g_hat_D_p_k = forwardResponseSensitivityComputer.create_D_g_hat_D_p();
778 }
779
780 // D.b) Zero out d_hat and D_d_hat_D_p if we are doing a summation type of
781 // evaluation
782 if (responseType_ == FSIAMET::RESPONSE_TYPE_SUM) {
783 if (!is_null(d_hat)) {
784 assign( d_hat.ptr(), ST::zero() );
785 }
786 if (!D_d_hat_D_p.isEmpty()) {
787 assign( D_d_hat_D_p.getMultiVector().ptr(), ST::zero() );
788 }
789 }
790
791 // D.c) Get product vector and multi-vector interfaces if
792 // we are just blocking up response functions
793 RCP<Thyra::ProductVectorBase<Scalar> > prod_d_hat;
794 RCP<Thyra::ProductMultiVectorBase<Scalar> > prod_D_d_hat_D_p_trans;
795 if (responseType_ == FSIAMET::RESPONSE_TYPE_BLOCK) {
796 prod_d_hat = rcp_dynamic_cast<Thyra::ProductVectorBase<Scalar> >(
797 d_hat, true);
798 prod_D_d_hat_D_p_trans = rcp_dynamic_cast<Thyra::ProductMultiVectorBase<Scalar> >(
799 D_d_hat_D_p.getMultiVector(), true);
800 }
801
802 //
803 // E) Run the integrator at the response time points and assemble
804 // the response function and/or the response function
805 // derivative
806 //
807
808 if (trace) *out << "\nStarting integration and assembly loop ...\n";
809
810 const int numResponseFunctions = responseTimes_.size();
811
812 for (int k = 0; k < numResponseFunctions; ++k ) {
813
814 OSTab os_tab(out);
815
816 const Scalar t = responseTimes_[k];
817
818 //
819 // E.1) Integrate up to t and get x_bar and x_bar_dot
820 //
821 // Note, x_bar and x_bar_dot may be the state or the state+sens!
822
823 if (trace)
824 *out << "\nIntegrating state (or state+sens) for response k = "
825 << k << " for t = " << t << " ...\n";
826
827 RCP<const Thyra::VectorBase<Scalar> > x_bar, x_bar_dot;
828
829 {
830 RYTHMOS_FUNC_TIME_MONITOR(
831 "Rythmos:ForwardSensitivityIntegratorAsModelEvaluator::evalModel: integrate"
832 );
833 get_fwd_x_and_x_dot( *integrator, t, &x_bar, &x_bar_dot );
834 }
835
836 if (print_norms) {
837 *out << "\n||x_bar||inf = " << norms_inf(*x_bar) << std::endl;
838 *out << "\n||x_bar_dot||inf = " << norms_inf(*x_bar_dot) << std::endl;
839 }
840 if (print_x) {
841 *out << "\nx_bar = " << *x_bar << std::endl;
842 *out << "\nx_bar_dot = " << *x_bar_dot << std::endl;
843 }
844
845 // Extra underlying state and sensitivities
846 RCP<const Thyra::VectorBase<Scalar> > x, x_dot;
847 RCP<const Thyra::MultiVectorBase<Scalar> > S, S_dot;
848 if (integrateStateAndSens) {
849 extractStateAndSens( x_bar, x_bar_dot, &x, &S, &x_dot, &S_dot );
850 }
851 else {
852 x = x_bar;
853 x_dot = x_bar_dot;
854 }
855
856 //
857 // E.2) Evaluate the response function g_k and its derivatives at this
858 // time point
859 //
860
861 if (trace)
862 *out << "\nEvaluating response function k = " << k << " at time t = " << t << " ...\n";
863
864 RCP<const Thyra::ModelEvaluator<Scalar> > responseFunc = responseFuncs_[k];
865
866 MEB::InArgs<Scalar> responseInArgs = responseFunc->createInArgs();
867 responseInArgs.setArgs(responseFuncInArgs_[k]);
868 responseInArgs.set_p(p_index_,p);
869
870 forwardResponseSensitivityComputer.resetResponseFunction(
871 responseFunc, responseInArgs
872 );
873
874 if (!is_null(D_g_hat_D_p_k)) {
875 forwardResponseSensitivityComputer.computeResponseAndSensitivity(
876 x_dot.get(), S_dot.get(), *x, *S, t, g_k.get(), &*D_g_hat_D_p_k
877 );
878 }
879 else {
880 forwardResponseSensitivityComputer.computeResponse(
881 x_dot.get(), *x, t, g_k.get()
882 );
883 }
884
885 //
886 // E.3) Assemble the final response functions and derivatives
887 //
888
889 // E.3.a) Assemble d_hat from g_k
890 if (!is_null(d_hat)) {
891 if (responseType_ == FSIAMET::RESPONSE_TYPE_SUM) {
892 if (trace) *out << "\nd_hat += g_k ...\n";
893 Vp_V( d_hat.ptr(), *g_k );
894 }
895 else if (responseType_ == FSIAMET::RESPONSE_TYPE_BLOCK) {
896 if (trace) *out << "\nd_hat["<<k<<"] = g_k ...\n";
897 assign( prod_d_hat->getNonconstVectorBlock(k).ptr(), *g_k );
898 }
899 }
900
901 // E.3.b) Assemble D_d_hat_Dp from D_g_hat_D_p_k
902 if (!D_d_hat_D_p.isEmpty()) {
903 if (responseType_ == FSIAMET::RESPONSE_TYPE_SUM) {
904 if (trace) *out << "\nD_d_hat_D_p += D_g_hat_D_p_k ...\n";
905 Vp_V( D_d_hat_D_p.getMultiVector().ptr(), *D_g_hat_D_p_k );
906 if (dumpSensitivities_ || print_x)
907 *out << "\nD_d_hat_D_p = "
908 << Teuchos::describe(*D_d_hat_D_p.getMultiVector(),Teuchos::VERB_EXTREME);
909 }
910 else if (responseType_ == FSIAMET::RESPONSE_TYPE_BLOCK) {
911 if (trace) *out << "\nD_d_hat_D_p["<<k<<"] = D_g_hat_D_p_k ...\n";
912 assign(
913 prod_D_d_hat_D_p_trans->getNonconstMultiVectorBlock(k).ptr(),
914 *D_g_hat_D_p_k
915 );
916 }
917 }
918
919 }
920
921 THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_END();
922
923}
924
925
926} // namespace Rythmos
927
928
929#endif // RYTHMOS_FORWARD_SENSITIVITY_INTEGRATOR_AS_MODEL_EVALUATOR_HPP
Concrete Thyra::ModelEvaluator subclass that turns a forward ODE/DAE with an observation into a param...
void initialize(const RCP< StepperBase< Scalar > > &stateStepper, const RCP< IntegratorBase< Scalar > > &stateIntegrator, const RCP< ForwardSensitivityStepper< Scalar > > &stateAndSensStepper, const RCP< IntegratorBase< Scalar > > &stateAndSensIntegrator, const Thyra::ModelEvaluatorBase::InArgs< Scalar > &stateAndSensInitCond, const Array< Scalar > &responseTimes, const Array< RCP< const Thyra::ModelEvaluator< Scalar > > > &responseFuncs, const Array< Thyra::ModelEvaluatorBase::InArgs< Scalar > > &responseFuncBasePoints, const ForwardSensitivityIntegratorAsModelEvaluatorTypes::EResponseType responseType)
Initalized with all objects needed to define evaluation.