Rythmos - Transient Integration for Differential Equations  Version of the Day
 All Classes Functions Variables Typedefs Pages
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 
52 namespace Rythmos {
53 
54 
301 template<class Scalar>
303  : virtual public ForwardSensitivityModelEvaluatorBase<Scalar>
304 {
305 public:
306 
309 
312 
328  void initializeStructure(
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 
359  void setStateIntegrator(
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 
514 private:
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 
529 private:
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 
584 template<class Scalar>
585 RCP<ForwardSensitivityImplicitModelEvaluator<Scalar> > forwardSensitivityImplicitModelEvaluator()
586 {
588 }
589 
590 template<class Scalar>
592  : p_index_(0), np_(-1)
593 {}
594 
595 
596 template<class Scalar>
597 RCP<const Thyra::ModelEvaluator<Scalar> >
599 {
600  return stateModel_;
601 }
602 
603 
604 template<class Scalar>
605 RCP<Thyra::ModelEvaluator<Scalar> >
607 {
608  return Teuchos::null;
609 }
610 
611 
612 template<class Scalar>
614 {
615  return p_index_;
616 }
617 
618 
619 template<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 
629 template<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 
715 template<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 
775 template<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 
785 template<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 
795 template<class Scalar>
796 RCP<const Thyra::DefaultMultiVectorProductVectorSpace<Scalar> >
798 {
799  return s_bar_space_;
800 }
801 
802 
803 template<class Scalar>
804 RCP<const Thyra::VectorSpaceBase<Scalar> >
806 {
807  return p_space_;
808 }
809 
810 
811 // Public functions overridden from ModelEvaulator
812 
813 
814 template<class Scalar>
815 RCP<const Thyra::VectorSpaceBase<Scalar> >
817 {
818  return s_bar_space_;
819 }
820 
821 
822 template<class Scalar>
823 RCP<const Thyra::VectorSpaceBase<Scalar> >
825 {
826  return f_sens_space_;
827 }
828 
829 
830 template<class Scalar>
831 Thyra::ModelEvaluatorBase::InArgs<Scalar>
833 {
834  return nominalValues_;
835 }
836 
837 
838 template<class Scalar>
839 RCP<Thyra::LinearOpWithSolveBase<Scalar> >
841 {
842  return Thyra::multiVectorLinearOpWithSolve<Scalar>();
843 }
844 
845 
846 template<class Scalar>
847 Thyra::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 
869 template<class Scalar>
870 Thyra::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 
893 template<class Scalar>
894 void 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 
1012 template<class Scalar>
1013 void 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 
1086 template<class Scalar>
1087 void 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 
1106 template<class Scalar>
1107 void 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 initializeStructureInitCondOnly(const RCP< const Thyra::ModelEvaluator< Scalar > > &stateModel, const RCP< const Thyra::VectorSpaceBase< Scalar > > &p_space)
virtual Teuchos::RCP< Thyra::NonlinearSolverBase< Scalar > > getNonconstSolver()=0
RCP< const Thyra::DefaultMultiVectorProductVectorSpace< Scalar > > get_s_bar_space() const
void initializeStructure(const RCP< const Thyra::ModelEvaluator< Scalar > > &stateModel, const int p_index)
Intialize the with the model structure.
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.
Base class for defining stepper functionality.
Abstract interface for time integrators.
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.
Base class mix-in interface for single-residual model evaluators.
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.
Mix-in interface all implicit stepper objects that accept a nonlinear solver to be used to compute th...
Forward sensitivity transient ModelEvaluator node interface class.
virtual Scalar get_coeff_x_dot() const =0