Rythmos - Transient Integration for Differential Equations Version of the Day
Loading...
Searching...
No Matches
Rythmos_ForwardSensitivityImplicitModelEvaluator.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_IMPLICIT_MODEL_EVALUATOR_HPP
30#define RYTHMOS_FORWARD_SENSITIVITY_IMPLICIT_MODEL_EVALUATOR_HPP
31
32
33#include "Rythmos_IntegratorBase.hpp"
34#include "Rythmos_ForwardSensitivityModelEvaluatorBase.hpp"
35#include "Rythmos_SolverAcceptingStepperBase.hpp"
36#include "Rythmos_SingleResidualModelEvaluator.hpp"
37#include "Thyra_ModelEvaluator.hpp" // Interface
38#include "Thyra_StateFuncModelEvaluatorBase.hpp" // Implementation
39#include "Thyra_DefaultProductVectorSpace.hpp"
40#include "Thyra_PhysicallyBlockedLinearOpWithSolveBase.hpp" // Interface
41#include "Thyra_DefaultBlockedTriangularLinearOpWithSolve.hpp" // Implementation
42#include "Thyra_ModelEvaluatorDelegatorBase.hpp"
43#include "Thyra_ModelEvaluatorHelpers.hpp"
44#include "Thyra_DefaultMultiVectorProductVectorSpace.hpp"
45#include "Thyra_DefaultMultiVectorProductVector.hpp"
46#include "Thyra_DefaultMultiVectorLinearOpWithSolve.hpp"
47#include "Thyra_MultiVectorStdOps.hpp"
48#include "Teuchos_implicit_cast.hpp"
49#include "Teuchos_Assert.hpp"
50
51
52namespace Rythmos {
53
54
301template<class Scalar>
303 : virtual public ForwardSensitivityModelEvaluatorBase<Scalar>
304{
305public:
306
309
312
329 const RCP<const Thyra::ModelEvaluator<Scalar> > &stateModel,
330 const int p_index
331 );
332
334 RCP<const Thyra::ModelEvaluator<Scalar> >
335 getStateModel() const;
336
338 RCP<Thyra::ModelEvaluator<Scalar> >
339 getNonconstStateModel() const;
340
342 int get_p_index() const;
343
360 const RCP<IntegratorBase<Scalar> > &stateIntegrator,
361 const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateBasePoint
362 );
363
404 Ptr<StepperBase<Scalar> > stateStepper,
405 bool forceUpToDateW
406 );
407
408// /** \brief Initialize full state for a single point in time.
409// *
410// * \param stateBasePoint [in] The base point
411// * <tt>(x_dot_tilde,x_tilde,t_tilde,...)</tt> for which the sensitivities
412// * will be computed and represented at time <tt>t_tilde</tt>. Any
413// * sensitivities that are needed must be computed at this same time point.
414// * The value of <tt>alpha</tt> and <tt>beta</tt> here are ignored.
415// *
416// * \param W_tilde [in,persisting] The composite state derivative matrix
417// * computed at the base point <tt>stateBasePoint</tt> with
418// * <tt>alpha=coeff_x_dot</tt> and <tt>beta=coeff_x</tt>. This argument can
419// * be null, in which case it can be computed internally if needed or not at
420// * all.
421// *
422// * \param coeff_x_dot [in] The value of <tt>alpha</tt> for which
423// * <tt>W_tilde</tt> was computed.
424// *
425// * \param coeff_x [in] The value of <tt>beta</tt> for which <tt>W_tilde</tt>
426// * was computed.
427// *
428// * \param DfDx_dot [in] The matrix <tt>d(f)/d(x_dot)</tt> computed at
429// * <tt>stateBasePoint</tt> if available. If this argument is null, then it
430// * will be computed internally if needed. The default value is
431// * <tt>Teuchos::null</tt>.
432// *
433// * \param DfDp [in] The matrix <tt>d(f)/d(p(p_index))</tt> computed at
434// * <tt>stateBasePoint</tt> if available. If this argument is null, then it
435// * will be computed internally if needed. The default value is
436// * <tt>Teuchos::null</tt>.
437// *
438// * <b>Preconditions:</b><ul>
439// * <li><tt>!is_null(W_tilde)</tt>
440// * </ul>
441// *
442// * This function must be called after <tt>intializeStructure()</tt> and
443// * before <tt>evalModel()</tt> is called. After this function is called,
444// * <tt>*this</tt> model is fully initialized and ready to compute the
445// * requested outputs.
446// */
447// void initializePointState(
448// const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateBasePoint,
449// const RCP<const Thyra::LinearOpWithSolveBase<Scalar> > &W_tilde,
450// const Scalar &coeff_x_dot,
451// const Scalar &coeff_x,
452// const RCP<const Thyra::LinearOpBase<Scalar> > &DfDx_dot = Teuchos::null,
453// const RCP<const Thyra::MultiVectorBase<Scalar> > &DfDp = Teuchos::null
454// );
455
456 // 2007/05/22: rabartl: ToDo: Add an InterpolationBufferBase
457 // stateInterpBuffer object to the initializeState(...) function that can be
458 // used to get x and x_dot at different points in time t. Then, modify the
459 // logic to recompute all of the needed matrices if t != t_base (as passed
460 // in through stateBasePoint). The values of x(t) and xdot(t) can then be
461 // gotten from the stateInterpBuffer object!
462
464
467
470 const RCP<const Thyra::ModelEvaluator<Scalar> >& stateModel,
471 const RCP<const Thyra::VectorSpaceBase<Scalar> >& p_space
472 );
473
475 RCP<const Thyra::DefaultMultiVectorProductVectorSpace<Scalar> >
476 get_s_bar_space() const;
477
479 RCP<const Thyra::VectorSpaceBase<Scalar> > get_p_sens_space() const;
480
482
485
487 RCP<const Thyra::VectorSpaceBase<Scalar> > get_x_space() const;
489 RCP<const Thyra::VectorSpaceBase<Scalar> > get_f_space() const;
491 Thyra::ModelEvaluatorBase::InArgs<Scalar> getNominalValues() const;
493 RCP<Thyra::LinearOpWithSolveBase<Scalar> > create_W() const;
495 Thyra::ModelEvaluatorBase::InArgs<Scalar> createInArgs() const;
496
498
501
503 void initializeState(
504 const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateBasePoint,
505 const RCP<const Thyra::LinearOpWithSolveBase<Scalar> > &W_tilde,
506 const Scalar &coeff_x_dot,
507 const Scalar &coeff_x,
508 const RCP<const Thyra::LinearOpBase<Scalar> > &DfDx_dot = Teuchos::null,
509 const RCP<const Thyra::MultiVectorBase<Scalar> > &DfDp = Teuchos::null
510 );
511
513
514private:
515
518
520 Thyra::ModelEvaluatorBase::OutArgs<Scalar> createOutArgsImpl() const;
522 void evalModelImpl(
523 const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs,
524 const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs
525 ) const;
526
528
529private:
530
531 // /////////////////////////
532 // Private data members
533
534 RCP<const Thyra::ModelEvaluator<Scalar> > stateModel_;
535 int p_index_;
536 RCP<const Thyra::VectorSpaceBase<Scalar> > p_space_;
537 int np_;
538
539 RCP<IntegratorBase<Scalar> > stateIntegrator_;
540
541 RCP<const Thyra::DefaultMultiVectorProductVectorSpace<Scalar> > s_bar_space_;
542 RCP<const Thyra::DefaultMultiVectorProductVectorSpace<Scalar> > f_sens_space_;
543
544 Thyra::ModelEvaluatorBase::InArgs<Scalar> nominalValues_;
545
546 mutable Thyra::ModelEvaluatorBase::InArgs<Scalar> stateBasePoint_;
547
548 mutable RCP<const Thyra::LinearOpWithSolveBase<Scalar> > W_tilde_;
549 mutable Scalar coeff_x_dot_;
550 mutable Scalar coeff_x_;
551 mutable RCP<const Thyra::LinearOpBase<Scalar> > DfDx_dot_;
552 mutable RCP<const Thyra::MultiVectorBase<Scalar> > DfDp_;
553
554 mutable RCP<Thyra::LinearOpWithSolveBase<Scalar> > W_tilde_compute_;
555 mutable RCP<Thyra::LinearOpBase<Scalar> > DfDx_dot_compute_;
556 mutable RCP<Thyra::MultiVectorBase<Scalar> > DfDp_compute_;
557
558 // /////////////////////////
559 // Private member functions
560
561 bool hasStateFuncParams() const { return p_index_ >= 0; }
562
563 void initializeStructureCommon(
564 const RCP<const Thyra::ModelEvaluator<Scalar> > &stateModel,
565 const int p_index,
566 const RCP<const Thyra::VectorSpaceBase<Scalar> > &p_space
567 );
568
569 void wrapNominalValuesAndBounds();
570
571 void computeDerivativeMatrices(
572 const Thyra::ModelEvaluatorBase::InArgs<Scalar> &point
573 ) const;
574
575};
576
577
578// /////////////////////////////////
579// Implementations
580
581
582// Constructors/Intializers/Accessors
583
584template<class Scalar>
585RCP<ForwardSensitivityImplicitModelEvaluator<Scalar> > forwardSensitivityImplicitModelEvaluator()
586{
588}
589
590template<class Scalar>
592 : p_index_(0), np_(-1)
593{}
594
595
596template<class Scalar>
597RCP<const Thyra::ModelEvaluator<Scalar> >
599{
600 return stateModel_;
601}
602
603
604template<class Scalar>
605RCP<Thyra::ModelEvaluator<Scalar> >
607{
608 return Teuchos::null;
609}
610
611
612template<class Scalar>
614{
615 return p_index_;
616}
617
618
619template<class Scalar>
621 const RCP<IntegratorBase<Scalar> > &stateIntegrator,
622 const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateBasePoint
623 )
624{
625 stateIntegrator_ = stateIntegrator.assert_not_null();
626 stateBasePoint_ = stateBasePoint;
627}
628
629template<class Scalar>
631 Ptr<StepperBase<Scalar> > stateStepper,
632 bool forceUpToDateW
633 )
634{
635 TEUCHOS_ASSERT( Teuchos::nonnull(stateStepper) );
636 Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
637 Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
638
639 const StepStatus<Scalar> stateStepStatus = stateStepper->getStepStatus();
640 TEUCHOS_TEST_FOR_EXCEPTION(
641 stateStepStatus.stepStatus != STEP_STATUS_CONVERGED, std::logic_error,
642 "Error, the status should be converged since a positive step size was returned!"
643 );
644 Scalar curr_t = stateStepStatus.time;
645
646 // Get both x and x_dot since these are needed compute other derivative
647 // objects at these points.
648 RCP<const Thyra::VectorBase<Scalar> > x, x_dot;
649 get_x_and_x_dot(*stateStepper,curr_t,&x,&x_dot);
650
651 stateBasePoint_ = stateStepper->getInitialCondition();
652 stateBasePoint_.set_x_dot( x_dot );
653 stateBasePoint_.set_x( x );
654 stateBasePoint_.set_t( curr_t );
655
656 // Grab the SingleResidualModel that was used to compute the state timestep.
657 // From this, we can get the constants that where used to compute W!
658 RCP<SolverAcceptingStepperBase<Scalar> >
659 sasStepper = Teuchos::rcp_dynamic_cast<SolverAcceptingStepperBase<Scalar> >(
660 Teuchos::rcpFromRef(*stateStepper),true
661 );
662 RCP<Thyra::NonlinearSolverBase<Scalar> >
663 stateTimeStepSolver = sasStepper->getNonconstSolver();
664 RCP<const SingleResidualModelEvaluatorBase<Scalar> >
665 singleResidualModel
666 = Teuchos::rcp_dynamic_cast<const SingleResidualModelEvaluatorBase<Scalar> >(
667 stateTimeStepSolver->getModel(), true
668 );
669 const Scalar
670 coeff_x_dot = singleResidualModel->get_coeff_x_dot(),
671 coeff_x = singleResidualModel->get_coeff_x();
672
673 // Get W (and force an update if not up to date already)
674
675 using Teuchos::as;
676 if ((as<int>(verbLevel) >= as<int>(Teuchos::VERB_MEDIUM)) && forceUpToDateW) {
677 *out << "\nForcing an update of W at the converged state timestep ...\n";
678 }
679
680 RCP<Thyra::LinearOpWithSolveBase<Scalar> >
681 W_tilde = stateTimeStepSolver->get_nonconst_W(forceUpToDateW);
682
683 TEUCHOS_TEST_FOR_EXCEPTION(
684 is_null(W_tilde), std::logic_error,
685 "Error, the W from the state time step must be non-null!"
686 );
687
688#ifdef HAVE_RYTHMOS_DEBUG
689 TEUCHOS_TEST_FOR_EXCEPTION(
690 is_null(stateModel_), std::logic_error,
691 "Error, you must call intializeStructure(...) before you call initializeState(...)"
692 );
693 TEUCHOS_TEST_FOR_EXCEPT( is_null(stateBasePoint_.get_x()) );
694 TEUCHOS_TEST_FOR_EXCEPT( is_null(stateBasePoint_.get_x_dot()) );
695 TEUCHOS_TEST_FOR_EXCEPT( is_null(stateBasePoint_.get_p(p_index_)) );
696 // What about the other parameter values? We really can't say anything
697 // about these and we can't check them. They can be null just fine.
698 if (!is_null(W_tilde)) {
699 THYRA_ASSERT_VEC_SPACES("",*W_tilde->domain(),*stateModel_->get_x_space());
700 THYRA_ASSERT_VEC_SPACES("",*W_tilde->range(),*stateModel_->get_f_space());
701 }
702#endif
703
704 W_tilde_ = W_tilde;
705 coeff_x_dot_ = coeff_x_dot;
706 coeff_x_ = coeff_x;
707 DfDx_dot_ = Teuchos::null;
708 DfDp_ = Teuchos::null;
709
710 wrapNominalValuesAndBounds();
711
712}
713
714
715template<class Scalar>
717 const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateBasePoint,
718 const RCP<const Thyra::LinearOpWithSolveBase<Scalar> > &W_tilde,
719 const Scalar &coeff_x_dot,
720 const Scalar &coeff_x,
721 const RCP<const Thyra::LinearOpBase<Scalar> > &DfDx_dot,
722 const RCP<const Thyra::MultiVectorBase<Scalar> > &DfDp
723 )
724{
725
726 // typedef Thyra::ModelEvaluatorBase MEB; // unused
727
728#ifdef HAVE_RYTHMOS_DEBUG
729 TEUCHOS_TEST_FOR_EXCEPTION(
730 is_null(stateModel_), std::logic_error,
731 "Error, you must call intializeStructure(...) before you call initializeState(...)"
732 );
733 TEUCHOS_TEST_FOR_EXCEPT( is_null(stateBasePoint.get_x()) );
734 TEUCHOS_TEST_FOR_EXCEPT( is_null(stateBasePoint.get_x_dot()) );
735 if (hasStateFuncParams()) {
736 TEUCHOS_TEST_FOR_EXCEPT( is_null(stateBasePoint.get_p(p_index_)) );
737 }
738 // What about the other parameter values? We really can't say anything
739 // about these and we can't check them. They can be null just fine.
740 if (nonnull(W_tilde)) {
741 THYRA_ASSERT_VEC_SPACES("", *W_tilde->domain(), *stateModel_->get_x_space());
742 THYRA_ASSERT_VEC_SPACES("", *W_tilde->range(), *stateModel_->get_f_space());
743 }
744 if (nonnull(DfDx_dot)) {
745 THYRA_ASSERT_VEC_SPACES("", *DfDx_dot->domain(), *stateModel_->get_x_space());
746 THYRA_ASSERT_VEC_SPACES("", *DfDx_dot->range(), *stateModel_->get_f_space());
747 }
748 if (nonnull(DfDp)) {
749 TEUCHOS_ASSERT(hasStateFuncParams());
750 THYRA_ASSERT_VEC_SPACES("", *DfDp->domain(), *p_space_);
751 THYRA_ASSERT_VEC_SPACES("", *DfDp->range(), *stateModel_->get_f_space());
752 }
753#endif
754
755 stateBasePoint_ = stateBasePoint;
756
757 // Set whatever derivatives where passed in. If an input in null, then the
758 // member will be null and the null linear operators will be computed later
759 // just in time.
760
761 W_tilde_ = W_tilde;
762 coeff_x_dot_ = coeff_x_dot;
763 coeff_x_ = coeff_x;
764 DfDx_dot_ = DfDx_dot;
765 DfDp_ = DfDp;
766
767 wrapNominalValuesAndBounds();
768
769}
770
771
772// Public functions overridden from ForwardSensitivityModelEvaluatorBase
773
774
775template<class Scalar>
777 const RCP<const Thyra::ModelEvaluator<Scalar> >& stateModel,
778 const int p_index
779 )
780{
781 initializeStructureCommon(stateModel, p_index, Teuchos::null);
782}
783
784
785template<class Scalar>
787 const RCP<const Thyra::ModelEvaluator<Scalar> >& stateModel,
788 const RCP<const Thyra::VectorSpaceBase<Scalar> >& p_space
789 )
790{
791 initializeStructureCommon(stateModel, -1, p_space);
792}
793
794
795template<class Scalar>
796RCP<const Thyra::DefaultMultiVectorProductVectorSpace<Scalar> >
798{
799 return s_bar_space_;
800}
801
802
803template<class Scalar>
804RCP<const Thyra::VectorSpaceBase<Scalar> >
806{
807 return p_space_;
808}
809
810
811// Public functions overridden from ModelEvaulator
812
813
814template<class Scalar>
815RCP<const Thyra::VectorSpaceBase<Scalar> >
817{
818 return s_bar_space_;
819}
820
821
822template<class Scalar>
823RCP<const Thyra::VectorSpaceBase<Scalar> >
825{
826 return f_sens_space_;
827}
828
829
830template<class Scalar>
831Thyra::ModelEvaluatorBase::InArgs<Scalar>
833{
834 return nominalValues_;
835}
836
837
838template<class Scalar>
839RCP<Thyra::LinearOpWithSolveBase<Scalar> >
841{
842 return Thyra::multiVectorLinearOpWithSolve<Scalar>();
843}
844
845
846template<class Scalar>
847Thyra::ModelEvaluatorBase::InArgs<Scalar>
849{
850 typedef Thyra::ModelEvaluatorBase MEB;
851 MEB::InArgs<Scalar> stateModelInArgs = stateModel_->createInArgs();
852 MEB::InArgsSetup<Scalar> inArgs;
853 inArgs.setModelEvalDescription(this->description());
854 inArgs.setSupports( MEB::IN_ARG_x_dot,
855 stateModelInArgs.supports(MEB::IN_ARG_x_dot) );
856 inArgs.setSupports( MEB::IN_ARG_x );
857 inArgs.setSupports( MEB::IN_ARG_t );
858 inArgs.setSupports( MEB::IN_ARG_alpha,
859 stateModelInArgs.supports(MEB::IN_ARG_alpha) );
860 inArgs.setSupports( MEB::IN_ARG_beta,
861 stateModelInArgs.supports(MEB::IN_ARG_beta) );
862 return inArgs;
863}
864
865
866// Private functions overridden from ModelEvaulatorDefaultBase
867
868
869template<class Scalar>
870Thyra::ModelEvaluatorBase::OutArgs<Scalar>
872{
873
874 typedef Thyra::ModelEvaluatorBase MEB;
875
876 MEB::OutArgs<Scalar> stateModelOutArgs = stateModel_->createOutArgs();
877 MEB::OutArgsSetup<Scalar> outArgs;
878
879 outArgs.setModelEvalDescription(this->description());
880
881 outArgs.setSupports(MEB::OUT_ARG_f);
882
883 if (stateModelOutArgs.supports(MEB::OUT_ARG_W) ) {
884 outArgs.setSupports(MEB::OUT_ARG_W);
885 outArgs.set_W_properties(stateModelOutArgs.get_W_properties());
886 }
887
888 return outArgs;
889
890}
891
892
893template<class Scalar>
894void ForwardSensitivityImplicitModelEvaluator<Scalar>::evalModelImpl(
895 const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs,
896 const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs
897 ) const
898{
899
900 using Teuchos::as;
901 using Teuchos::rcp_dynamic_cast;
902 typedef Teuchos::ScalarTraits<Scalar> ST;
903 // typedef Thyra::ModelEvaluatorBase MEB; // unused
904 typedef Teuchos::VerboseObjectTempState<Thyra::ModelEvaluatorBase> VOTSME;
905
906 THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_GEN_BEGIN(
907 "ForwardSensitivityImplicitModelEvaluator", inArgs, outArgs, Teuchos::null );
908
909 //
910 // Update the derivative matrices if they are not already updated for the
911 // given time!.
912 //
913
914 {
915 RYTHMOS_FUNC_TIME_MONITOR_DIFF(
916 "Rythmos:ForwardSensitivityImplicitModelEvaluator::evalModel: computeMatrices",
917 RythmosFSIMEmain
918 );
919 computeDerivativeMatrices(inArgs);
920 }
921
922 //
923 // InArgs
924 //
925
926 RCP<const Thyra::DefaultMultiVectorProductVector<Scalar> >
927 s_bar = rcp_dynamic_cast<const Thyra::DefaultMultiVectorProductVector<Scalar> >(
928 inArgs.get_x().assert_not_null(), true
929 );
930 RCP<const Thyra::DefaultMultiVectorProductVector<Scalar> >
931 s_bar_dot = rcp_dynamic_cast<const Thyra::DefaultMultiVectorProductVector<Scalar> >(
932 inArgs.get_x_dot().assert_not_null(), true
933 );
934 const Scalar
935 alpha = inArgs.get_alpha();
936 const Scalar
937 beta = inArgs.get_beta();
938
939 RCP<const Thyra::MultiVectorBase<Scalar> >
940 S = s_bar->getMultiVector();
941 RCP<const Thyra::MultiVectorBase<Scalar> >
942 S_dot = s_bar_dot->getMultiVector();
943
944 //
945 // OutArgs
946 //
947
948 RCP<Thyra::DefaultMultiVectorProductVector<Scalar> >
949 f_sens = rcp_dynamic_cast<Thyra::DefaultMultiVectorProductVector<Scalar> >(
950 outArgs.get_f(), true
951 );
952 RCP<Thyra::DefaultMultiVectorLinearOpWithSolve<Scalar> >
953 W_sens = rcp_dynamic_cast<Thyra::DefaultMultiVectorLinearOpWithSolve<Scalar> >(
954 outArgs.get_W(), true
955 );
956
957 RCP<Thyra::MultiVectorBase<Scalar> >
958 F_sens = f_sens->getNonconstMultiVector().assert_not_null();
959
960 //
961 // Compute the requested functions
962 //
963
964 if(nonnull(F_sens)) {
965
966 RYTHMOS_FUNC_TIME_MONITOR_DIFF(
967 "Rythmos:ForwardSensitivityImplicitModelEvaluator::evalModel: computeSens",
968 Rythmos_FSIME);
969
970 // S_diff = -(coeff_x_dot/coeff_x)*S + S_dot
971 RCP<Thyra::MultiVectorBase<Scalar> >
972 S_diff = createMembers( stateModel_->get_x_space(), np_ );
973 V_StVpV( S_diff.ptr(), as<Scalar>(-coeff_x_dot_/coeff_x_), *S, *S_dot );
974 // F_sens = (1/coeff_x) * W_tilde * S
975 Thyra::apply(
976 *W_tilde_, Thyra::NOTRANS,
977 *S, F_sens.ptr(),
978 as<Scalar>(1.0/coeff_x_), ST::zero()
979 );
980 // F_sens += d(f)/d(x_dot) * S_diff
981 Thyra::apply(
982 *DfDx_dot_, Thyra::NOTRANS,
983 *S_diff, F_sens.ptr(),
984 ST::one(), ST::one()
985 );
986 // F_sens += d(f)/d(p)
987 if (hasStateFuncParams())
988 Vp_V( F_sens.ptr(), *DfDp_ );
989 }
990
991 if(nonnull(W_sens)) {
992 TEUCHOS_TEST_FOR_EXCEPTION(
993 alpha != coeff_x_dot_, std::logic_error,
994 "Error, alpha="<<alpha<<" != coeff_x_dot="<<coeff_x_dot_
995 <<" with difference = "<<(alpha-coeff_x_dot_)<<"!" );
996 TEUCHOS_TEST_FOR_EXCEPTION(
997 beta != coeff_x_, std::logic_error,
998 "Error, beta="<<beta<<" != coeff_x="<<coeff_x_
999 <<" with difference = "<<(beta-coeff_x_)<<"!" );
1000 W_sens->initialize( W_tilde_, s_bar_space_, f_sens_space_ );
1001
1002 }
1003
1004 THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_END();
1005
1006}
1007
1008
1009// private
1010
1011
1012template<class Scalar>
1013void ForwardSensitivityImplicitModelEvaluator<Scalar>::initializeStructureCommon(
1014 const RCP<const Thyra::ModelEvaluator<Scalar> >& stateModel,
1015 const int p_index,
1016 const RCP<const Thyra::VectorSpaceBase<Scalar> >& p_space
1017 )
1018{
1019
1020 typedef Thyra::ModelEvaluatorBase MEB;
1021
1022 //
1023 // Validate input
1024 //
1025
1026 TEUCHOS_ASSERT(nonnull(stateModel));
1027 TEUCHOS_ASSERT(p_index >= 0 || nonnull(p_space));
1028 if (p_index >= 0) {
1029 TEUCHOS_TEST_FOR_EXCEPTION(
1030 !( 0 <= p_index && p_index < stateModel->Np() ), std::logic_error,
1031 "Error, p_index does not fall in the range [0,"<<(stateModel->Np()-1)<<"]!" );
1032 // ToDo: Validate support for DfDp!
1033 }
1034 else {
1035 TEUCHOS_ASSERT_EQUALITY(p_index, -1);
1036 }
1037
1038 //
1039 // Set the input objects
1040 //
1041
1042 stateModel_ = stateModel;
1043
1044 if (p_index >= 0) {
1045 p_index_ = p_index;
1046 p_space_ = stateModel_->get_p_space(p_index);
1047 }
1048 else {
1049 p_index_ = -1;
1050 p_space_ = p_space;
1051 }
1052
1053 np_ = p_space_->dim();
1054
1055 //
1056 // Create the structure of the model
1057 //
1058
1059 s_bar_space_ = Thyra::multiVectorProductVectorSpace(
1060 stateModel_->get_x_space(), np_
1061 );
1062
1063 f_sens_space_ = Thyra::multiVectorProductVectorSpace(
1064 stateModel_->get_f_space(), np_
1065 );
1066
1067 nominalValues_ = this->createInArgs();
1068
1069 //
1070 // Wipe out matrix storage
1071 //
1072
1073 stateBasePoint_ = MEB::InArgs<Scalar>();
1074 W_tilde_ = Teuchos::null;
1075 coeff_x_dot_ = 0.0;
1076 coeff_x_ = 0.0;
1077 DfDx_dot_ = Teuchos::null;
1078 DfDp_ = Teuchos::null;
1079 W_tilde_compute_ = Teuchos::null;
1080 DfDx_dot_compute_ = Teuchos::null;
1081 DfDp_compute_ = Teuchos::null;
1082
1083}
1084
1085
1086template<class Scalar>
1087void ForwardSensitivityImplicitModelEvaluator<Scalar>::wrapNominalValuesAndBounds()
1088{
1089
1090 // using Teuchos::rcp_dynamic_cast; // unused
1091 // typedef Thyra::ModelEvaluatorBase MEB; // unused
1092
1093 // nominalValues_.clear(); // ToDo: Implement this!
1094
1095 nominalValues_.set_t(stateModel_->getNominalValues().get_t());
1096
1097 // 2007/05/22: rabartl: Note: Currently there is not much of a reason to set
1098 // an initial condition here since the initial condition for the
1099 // sensitivities is really being set in the ForwardSensitivityStepper
1100 // object! In the future, a more general use of this class might benefit
1101 // from setting the initial condition here.
1102
1103}
1104
1105
1106template<class Scalar>
1107void ForwardSensitivityImplicitModelEvaluator<Scalar>::computeDerivativeMatrices(
1108 const Thyra::ModelEvaluatorBase::InArgs<Scalar> &point
1109 ) const
1110{
1111
1112 typedef Thyra::ModelEvaluatorBase MEB;
1113 typedef Teuchos::VerboseObjectTempState<MEB> VOTSME;
1114
1115 Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
1116 Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
1117
1118 const Scalar
1119 t_base = stateBasePoint_.get_t(),
1120 t = point.get_t();
1121
1122 TEUCHOS_ASSERT_EQUALITY( t , t_base );
1123
1124 if (is_null(W_tilde_)) {
1125 TEUCHOS_TEST_FOR_EXCEPT_MSG(true, "ToDo: compute W_tilde from scratch!");
1126 }
1127
1128 if ( is_null(DfDx_dot_) || is_null(DfDp_) ) {
1129
1130 MEB::InArgs<Scalar> inArgs = stateBasePoint_;
1131 MEB::OutArgs<Scalar> outArgs = stateModel_->createOutArgs();
1132
1133 Teuchos::RCP<Thyra::LinearOpBase<Scalar> > DfDx_dot_compute;
1134 if (is_null(DfDx_dot_)) {
1135 DfDx_dot_compute = stateModel_->create_W_op();
1136 inArgs.set_alpha(1.0);
1137 inArgs.set_beta(0.0);
1138 outArgs.set_W_op(DfDx_dot_compute);
1139 }
1140
1141 Teuchos::RCP<Thyra::MultiVectorBase<Scalar> > DfDp_compute;
1142 if (is_null(DfDp_) && hasStateFuncParams()) {
1143 DfDp_compute = Thyra::create_DfDp_mv(
1144 *stateModel_, p_index_,
1145 MEB::DERIV_MV_BY_COL
1146 ).getMultiVector();
1147 outArgs.set_DfDp(
1148 p_index_,
1149 MEB::Derivative<Scalar>(DfDp_compute,MEB::DERIV_MV_BY_COL)
1150 );
1151 }
1152
1153 VOTSME stateModel_outputTempState(stateModel_,out,verbLevel);
1154 stateModel_->evalModel(inArgs,outArgs);
1155 if (nonnull(DfDx_dot_compute))
1156 DfDx_dot_ = DfDx_dot_compute;
1157 if (nonnull(DfDp_compute))
1158 DfDp_ = DfDp_compute;
1159
1160 }
1161
1162 TEUCHOS_TEST_FOR_EXCEPT_MSG( nonnull(stateIntegrator_),
1163 "ToDo: Update for using the stateIntegrator!" );
1164
1165/* 2007/12/11: rabartl: ToDo: Update the code below to work for the general
1166 * case. It does not work in its current form!
1167 */
1168
1169/*
1170
1171 typedef Thyra::ModelEvaluatorBase MEB;
1172 typedef Teuchos::VerboseObjectTempState<MEB> VOTSME;
1173
1174 RCP<Teuchos::FancyOStream> out = this->getOStream();
1175 Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
1176
1177 const Scalar t = point.get_t();
1178
1179 //
1180 // A) Set up the base point at time t and determine what needs to be
1181 // computed here.
1182 //
1183
1184 bool update_W_tilde = false;
1185 bool update_DfDx_dot = false;
1186 bool update_DfDp = false;
1187
1188 if (nonnull(stateIntegrator_)) {
1189 // Get x and x_dot at t and flag to recompute all matrices!
1190 RCP<const Thyra::VectorBase<Scalar> > x, x_dot;
1191 get_fwd_x_and_x_dot( *stateIntegrator_, t, &x, &x_dot );
1192 stateBasePoint_.set_t(t);
1193 stateBasePoint_.set_x(x);
1194 stateBasePoint_.set_x_dot(x_dot);
1195 update_W_tilde = true;
1196 update_DfDx_dot = true;
1197 update_DfDp = true;
1198 }
1199 else {
1200 // The time point should be the same so only compute matrices that have
1201 // not already been computed.
1202 TEUCHOS_ASSERT_EQUALITY( t , stateBasePoint_.get_t() );
1203 if (is_null(W_tilde_))
1204 update_W_tilde = true;
1205 if (is_null(DfDx_dot_))
1206 update_DfDx_dot = true;
1207 if (is_null(DfDp_))
1208 update_DfDx_dot = true;
1209 }
1210
1211 //
1212 // B) Compute DfDx_dot and DfDp at the same time if needed
1213 //
1214
1215 if ( update_DfDx_dot || update_DfDp ) {
1216
1217 // B.1) Allocate storage for matrices. Note: All of these matrices are
1218 // needed so any of these that is null must be coputed!
1219
1220 if ( is_null(DfDx_dot_) || is_null(DfDx_dot_compute_) )
1221 DfDx_dot_compute_ = stateModel_->create_W_op();
1222
1223 if ( is_null(DfDp_) || is_null(DfDp_compute_) )
1224 DfDp_compute_ = Thyra::create_DfDp_mv(
1225 *stateModel_,p_index_,
1226 MEB::DERIV_MV_BY_COL
1227 ).getMultiVector();
1228
1229 // B.2) Setup the InArgs and OutArgs
1230
1231 MEB::InArgs<Scalar> inArgs = stateBasePoint_;
1232 MEB::OutArgs<Scalar> outArgs = stateModel_->createOutArgs();
1233
1234 if (update_DfDx_dot) {
1235 inArgs.set_alpha(1.0);
1236 inArgs.set_beta(0.0);
1237 outArgs.set_W_op(DfDx_dot_compute_);
1238 }
1239 // 2007/12/07: rabartl: ToDo: I need to change the structure of the
1240 // derivative properties in OutArgs to keep track of whether DfDx_dot is
1241 // constant or not separate from W. Right now I can't do this. This is
1242 // one reason to add a new DfDx_dot output object to OutArgs. A better
1243 // solution is a more general data structure giving the dependence of f()
1244 // and g[j]() on x, x_dot, p[l], t, etc ...
1245
1246 if (update_DfDp) {
1247 outArgs.set_DfDp(
1248 p_index_,
1249 MEB::Derivative<Scalar>(DfDp_compute_, MEB::DERIV_MV_BY_COL)
1250 );
1251 }
1252
1253 // B.3) Evaluate the outputs
1254
1255 VOTSME stateModel_outputTempState(stateModel_,out,verbLevel);
1256 stateModel_->evalModel(inArgs,outArgs);
1257
1258 // B.4) Set the outputs
1259
1260 if (nonnull(DfDx_dot_compute_))
1261 DfDx_dot_ = DfDx_dot_compute_;
1262
1263 if (nonnull(DfDp_compute_))
1264 DfDp_ = DfDp_compute_;
1265
1266 }
1267
1268 //
1269 // C) Compute W_tilde separately if needed
1270 //
1271
1272 if ( update_W_tilde ) {
1273
1274 // C.1) Allocate storage for matrices. Note: All of these matrices are
1275 // needed so any of these that is null must be coputed!
1276
1277 if ( is_null(W_tilde_) || is_null(W_tilde_compute_) )
1278 W_tilde_compute_ = stateModel_->create_W();
1279
1280 // C.2) Setup the InArgs and OutArgs
1281
1282 MEB::InArgs<Scalar> inArgs = stateBasePoint_;
1283 MEB::OutArgs<Scalar> outArgs = stateModel_->createOutArgs();
1284
1285 if (is_null(W_tilde_)) {
1286 coeff_x_dot_ = point.get_alpha();
1287 coeff_x_ = point.get_beta();
1288 inArgs.set_alpha(coeff_x_dot_);
1289 inArgs.set_beta(coeff_x_);
1290 outArgs.set_W(W_tilde_compute_);
1291 }
1292
1293 // C.3) Evaluate the outputs
1294
1295 VOTSME stateModel_outputTempState(stateModel_,out,verbLevel);
1296 stateModel_->evalModel(inArgs,outArgs);
1297
1298 // B.4) Set the outputs
1299
1300 if (nonnull(W_tilde_compute_))
1301 W_tilde_ = W_tilde_compute_;
1302
1303 }
1304
1305 // 2007/12/07: rabartl: Note: Above, we see an example of where it would be
1306 // good to have a separate OutArg for DfDx_dot (as a LOB object) so that we
1307 // can compute W and DfDx_dot (and any other output quantitiy) at the same
1308 // time.
1309
1310*/
1311
1312
1313}
1314
1315
1316} // namespace Rythmos
1317
1318
1319#endif // RYTHMOS_FORWARD_SENSITIVITY_IMPLICIT_MODEL_EVALUATOR_HPP
void initializeState(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &stateBasePoint, const RCP< const Thyra::LinearOpWithSolveBase< Scalar > > &W_tilde, const Scalar &coeff_x_dot, const Scalar &coeff_x, const RCP< const Thyra::LinearOpBase< Scalar > > &DfDx_dot=Teuchos::null, const RCP< const Thyra::MultiVectorBase< Scalar > > &DfDp=Teuchos::null)
Deprecated.
void initializeStructure(const RCP< const Thyra::ModelEvaluator< Scalar > > &stateModel, const int p_index)
Intialize the with the model structure.
void setStateIntegrator(const RCP< IntegratorBase< Scalar > > &stateIntegrator, const Thyra::ModelEvaluatorBase::InArgs< Scalar > &stateBasePoint)
Set the state integrator that will be used to get x and x_dot at various time points.
RCP< const Thyra::VectorSpaceBase< Scalar > > get_p_sens_space() const
void initializePointState(Ptr< StepperBase< Scalar > > stateStepper, bool forceUpToDateW)
Initialize full state for a single point in time.
void initializeStructureInitCondOnly(const RCP< const Thyra::ModelEvaluator< Scalar > > &stateModel, const RCP< const Thyra::VectorSpaceBase< Scalar > > &p_space)
RCP< const Thyra::DefaultMultiVectorProductVectorSpace< Scalar > > get_s_bar_space() const
Forward sensitivity transient ModelEvaluator node interface class.
Abstract interface for time integrators.
Base class for defining stepper functionality.