Rythmos - Transient Integration for Differential Equations  Version of the Day
 All Classes Functions Variables Typedefs Pages
Rythmos_ForwardSensitivityStepper.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_STEPPER_HPP
30 #define RYTHMOS_FORWARD_SENSITIVITY_STEPPER_HPP
31 
32 
33 #include "Rythmos_StepperBase.hpp"
34 #include "Rythmos_StepperHelpers.hpp"
35 #include "Rythmos_ForwardSensitivityModelEvaluatorBase.hpp"
36 #include "Rythmos_ForwardSensitivityImplicitModelEvaluator.hpp"
37 #include "Rythmos_ForwardSensitivityExplicitModelEvaluator.hpp"
38 #include "Rythmos_StateAndForwardSensitivityModelEvaluator.hpp"
39 #include "Rythmos_SolverAcceptingStepperBase.hpp"
40 #include "Rythmos_IntegratorBase.hpp"
41 #include "Rythmos_SingleResidualModelEvaluatorBase.hpp"
42 #include "Thyra_ModelEvaluatorHelpers.hpp"
43 #include "Thyra_LinearNonlinearSolver.hpp"
44 #include "Thyra_ProductVectorBase.hpp"
45 #include "Thyra_AssertOp.hpp"
46 #include "Teuchos_ParameterListAcceptorDefaultBase.hpp"
47 #include "Teuchos_VerboseObjectParameterListHelpers.hpp"
48 #include "Teuchos_ConstNonconstObjectContainer.hpp"
49 #include "Teuchos_Assert.hpp"
50 #include "Teuchos_as.hpp"
51 
52 
53 namespace Rythmos {
54 
55 
208 template<class Scalar>
210  : virtual public StepperBase<Scalar>,
211  virtual public Teuchos::ParameterListAcceptorDefaultBase
212 {
213 public:
214 
216  typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType ScalarMag;
217 
220 
223 
278  const RCP<const Thyra::ModelEvaluator<Scalar> > &stateModel,
279  const int p_index,
280  const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateBasePoint,
281  const RCP<StepperBase<Scalar> > &stateStepper,
282  const RCP<Thyra::NonlinearSolverBase<Scalar> > &stateTimeStepSolver,
283  const RCP<StepperBase<Scalar> > &sensStepper = Teuchos::null,
284  const RCP<Thyra::NonlinearSolverBase<Scalar> > &sensTimeStepSolver = Teuchos::null
285  );
286 
312  const RCP<const Thyra::ModelEvaluator<Scalar> >& stateModel,
313  const RCP<const Thyra::VectorSpaceBase<Scalar> >& p_space,
314  const Thyra::ModelEvaluatorBase::InArgs<Scalar>& stateBasePoint,
315  const RCP<StepperBase<Scalar> >& stateStepper,
316  const RCP<Thyra::NonlinearSolverBase<Scalar> >& stateTimeStepSolver,
317  const RCP<StepperBase<Scalar> >& sensStepper = Teuchos::null,
318  const RCP<Thyra::NonlinearSolverBase<Scalar> >& sensTimeStepSolver = Teuchos::null
319  );
320 
357  const RCP<const Thyra::ModelEvaluator<Scalar> > &stateModel,
358  const int p_index,
359  const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateBasePoint,
360  const RCP<StepperBase<Scalar> > &stateStepper,
361  const RCP<Thyra::NonlinearSolverBase<Scalar> > &stateTimeStepSolver,
362  const RCP<IntegratorBase<Scalar> > &stateIntegrator,
363  const Scalar &finalTime,
364  const RCP<StepperBase<Scalar> > &sensStepper = Teuchos::null,
365  const RCP<Thyra::NonlinearSolverBase<Scalar> > &sensTimeStepSolver = Teuchos::null
366  );
367 
369  bool stateModelIsConst() const;
370 
374  RCP<const Thyra::ModelEvaluator<Scalar> >
375  getStateModel() const;
376 
380  RCP<StepperBase<Scalar> >
382 
386  RCP<const ForwardSensitivityModelEvaluatorBase<Scalar> >
387  getFwdSensModel() const;
388 
395  RCP<const StateAndForwardSensitivityModelEvaluator<Scalar> >
396  getStateAndFwdSensModel() const;
397 
399 
402 
404  void setParameterList(RCP<Teuchos::ParameterList> const& paramList);
406  RCP<const Teuchos::ParameterList> getValidParameters() const;
407 
409 
412 
414  bool acceptsModel() const;
415 
417  void setModel(
418  const RCP<const Thyra::ModelEvaluator<Scalar> >& model
419  );
420 
422  void setNonconstModel(
423  const RCP<Thyra::ModelEvaluator<Scalar> >& model
424  );
425 
432  RCP<const Thyra::ModelEvaluator<Scalar> > getModel() const;
433 
435  RCP<Thyra::ModelEvaluator<Scalar> > getNonconstModel();
436 
437  // RAB: ToDo: 2007/05/15: I need to talk with Todd about potentially
438  // removing the setModel() and getModel() functions from the StepperBase
439  // interface. In the case of this forward sensitivity solver, I am not sure
440  // that it makes a lot of sense to define a model. This surely will not be
441  // the model that a generic client would expect. The assumption I am sure
442  // would be that this model has the same space for x as the interpolation
443  // buffer but that is not true in this case.
444 
458  void setInitialCondition(
459  const Thyra::ModelEvaluatorBase::InArgs<Scalar> &state_and_sens_ic
460  );
461 
463  Thyra::ModelEvaluatorBase::InArgs<Scalar> getInitialCondition() const;
464 
466  Scalar takeStep( Scalar dt, StepSizeType stepType );
467 
469  const StepStatus<Scalar> getStepStatus() const;
470 
472 
475 
482  RCP<const Thyra::VectorSpaceBase<Scalar> >
483  get_x_space() const;
484 
486  void addPoints(
487  const Array<Scalar>& time_vec,
488  const Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >& x_vec,
489  const Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >& xdot_vec
490  );
491 
494 
496  void getPoints(
497  const Array<Scalar>& time_vec,
498  Array<RCP<const Thyra::VectorBase<Scalar> > >* x_vec,
499  Array<RCP<const Thyra::VectorBase<Scalar> > >* xdot_vec,
500  Array<ScalarMag>* accuracy_vec
501  ) const;
502 
504  void getNodes(Array<Scalar>* time_vec) const;
505 
507  void removeNodes(Array<Scalar>& time_vec);
508 
510  int getOrder() const;
511 
513 
516 
519  const RCP<const Thyra::ModelEvaluator<Scalar> > &stateModel,
520  const int p_index,
521  const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateBasePoint,
522  const RCP<StepperBase<Scalar> > &stateStepper,
523  const RCP<Thyra::NonlinearSolverBase<Scalar> > &stateTimeStepSolver,
524  const RCP<StepperBase<Scalar> > &sensStepper = Teuchos::null,
525  const RCP<Thyra::NonlinearSolverBase<Scalar> > &sensTimeStepSolver = Teuchos::null
526  )
527  {
529  stateModel, p_index, stateBasePoint, stateStepper, stateTimeStepSolver,
530  sensStepper, sensTimeStepSolver
531  );
532  }
533 
535 
536 private:
537  // ///////////////////
538  // Private types
539 
540  typedef Teuchos::ConstNonconstObjectContainer<Thyra::ModelEvaluator<Scalar> > CNCME;
541 
542  // /////////////////////////
543  // Private data members
544 
545  bool forceUpToDateW_;
546  CNCME stateModel_;
547  Thyra::ModelEvaluatorBase::InArgs<Scalar> stateBasePoint_;
548  RCP<StepperBase<Scalar> > stateStepper_;
549  RCP<Thyra::NonlinearSolverBase<Scalar> > stateTimeStepSolver_;
550  RCP<IntegratorBase<Scalar> > stateIntegrator_;
551  Scalar finalTime_;
552  Thyra::ModelEvaluatorBase::InArgs<Scalar> stateAndSensBasePoint_;
553  RCP<StepperBase<Scalar> > sensStepper_;
554  RCP<Thyra::NonlinearSolverBase<Scalar> > sensTimeStepSolver_;
555 
556  bool isSingleResidualStepper_;
557  RCP<ForwardSensitivityModelEvaluatorBase<Scalar> > sensModel_;
558  RCP<StateAndForwardSensitivityModelEvaluator<Scalar> > stateAndSensModel_;
559  Thyra::ModelEvaluatorBase::InArgs<Scalar> stateBasePoint_t_;
560 
561  static const std::string forceUpToDateW_name_;
562  static const bool forceUpToDateW_default_;
563 
564  // /////////////////////////
565  // Private member functions
566 
567 
568  // Common initialization helper
569  //
570  // Preconditions:
571  // (*) p_index >=0 or nonnull(p_space) == true
572  //
573  void initializeCommon(
574  const RCP<const Thyra::ModelEvaluator<Scalar> > &stateModel,
575  const int p_index,
576  const RCP<const Thyra::VectorSpaceBase<Scalar> > &p_space,
577  const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateBasePoint,
578  const RCP<StepperBase<Scalar> > &stateStepper,
579  const RCP<Thyra::NonlinearSolverBase<Scalar> > &stateTimeStepSolver,
580  const RCP<StepperBase<Scalar> > &sensStepper,
581  const RCP<Thyra::NonlinearSolverBase<Scalar> > &sensTimeStepSolver
582  );
583 
584  Scalar takeSyncedStep( Scalar dt, StepSizeType stepType );
585 
586  Scalar takeDecoupledStep( Scalar dt, StepSizeType stepType );
587 
588 };
589 
590 
591 // 2009/09/05: rabartl: ToDo: To fix the const and non-const handling of the
592 // stateModel in this class is going to be a lot of work but here is what
593 // needs to be done:
594 //
595 // (*) Duplicate each function that sets the stateModel, one for const and one
596 // for non-const.
597 //
598 // (*) Create a single a private version for each of these functions that
599 // accepts a Teuchos::ConstNonconstObjectContainer<> object and will implement
600 // the guts of the set up same as the existing functions.
601 //
602 // (*) Get all of the concrete StepperBase subclasses to implement the
603 // setModel(const RCP<const ME>&) and modelIsConst() functions and get them to
604 // use the Teuchos::ConstNonconstObjectContainer<> class as described above.
605 // This should be pretty easy as the only function that needs to be addressed
606 // in most cases is just the setModel(...) function.
607 //
608 
609 
614 template<class Scalar>
615 inline
616 RCP<ForwardSensitivityStepper<Scalar> >
618 {
619  return Teuchos::rcp(new ForwardSensitivityStepper<Scalar>());
620 }
621 
622 
627 template<class Scalar>
628 inline
629 RCP<ForwardSensitivityStepper<Scalar> >
631  const RCP<const Thyra::ModelEvaluator<Scalar> > &stateModel,
632  const int p_index,
633  const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateBasePoint,
634  const RCP<StepperBase<Scalar> > &stateStepper,
635  const RCP<Thyra::NonlinearSolverBase<Scalar> > &stateTimeStepSolver,
636  const RCP<StepperBase<Scalar> > &sensStepper = Teuchos::null,
637  const RCP<Thyra::NonlinearSolverBase<Scalar> > &sensTimeStepSolver = Teuchos::null
638  )
639 {
640  RCP<ForwardSensitivityStepper<Scalar> >
641  fwdSensStepper = Teuchos::rcp(new ForwardSensitivityStepper<Scalar>());
642  fwdSensStepper->initializeSyncedSteppers(
643  stateModel, p_index, stateBasePoint, stateStepper, stateTimeStepSolver );
644  return fwdSensStepper;
645 }
646 
647 
653 template<class Scalar>
655  const ForwardSensitivityStepper<Scalar> &fwdSensStepper
656  )
657 {
658  return fwdSensStepper.getFwdSensModel()->get_p_index();
659 }
660 
661 
668 template<class Scalar>
669 inline
670 Thyra::ModelEvaluatorBase::InArgs<Scalar>
672  const ForwardSensitivityStepper<Scalar> &fwdSensStepper,
673  const Thyra::ModelEvaluatorBase::InArgs<Scalar> &state_ic,
674  const RCP<const Thyra::MultiVectorBase<Scalar> > S_init = Teuchos::null,
675  const RCP<const Thyra::MultiVectorBase<Scalar> > S_dot_init = Teuchos::null
676  )
677 {
678 
679  using Teuchos::outArg;
680  using Thyra::assign;
681  typedef Thyra::ModelEvaluatorBase MEB;
682 
683  RCP<const Thyra::VectorBase<Scalar> > s_bar_init;
684  if (nonnull(S_init)) {
685  s_bar_init = create_s_bar_given_S(*fwdSensStepper.getFwdSensModel(), S_init);
686  }
687  else {
688  RCP<Thyra::VectorBase<Scalar> > s_bar_init_loc =
689  createMember(fwdSensStepper.getFwdSensModel()->get_x_space());
690  assign( outArg(*s_bar_init_loc), 0.0 );
691  s_bar_init = s_bar_init_loc;
692  }
693 
694  RCP<const Thyra::VectorBase<Scalar> > s_bar_dot_init;
695  if (nonnull(S_dot_init)) {
696  s_bar_dot_init = create_s_bar_given_S(*fwdSensStepper.getFwdSensModel(), S_dot_init);
697  }
698  else {
699  RCP<Thyra::VectorBase<Scalar> > s_bar_dot_init_loc =
700  createMember(fwdSensStepper.getFwdSensModel()->get_x_space());
701  assign( outArg(*s_bar_dot_init_loc), 0.0 );
702  s_bar_dot_init = s_bar_dot_init_loc;
703  }
704 
705  RCP<const Rythmos::StateAndForwardSensitivityModelEvaluator<Scalar> >
706  stateAndSensModel = fwdSensStepper.getStateAndFwdSensModel();
707 
708  MEB::InArgs<Scalar>
709  state_and_sens_ic = fwdSensStepper.getModel()->createInArgs();
710 
711  // Copy time, parameters etc.
712  state_and_sens_ic.setArgs(state_ic);
713  // Set initial condition for x_bar = [ x; s_bar ]
714  state_and_sens_ic.set_x(
715  stateAndSensModel->create_x_bar_vec(state_ic.get_x(), s_bar_init)
716  );
717  // Set initial condition for x_bar_dot = [ x_dot; s_bar_dot ]
718  state_and_sens_ic.set_x_dot(
719  stateAndSensModel->create_x_bar_vec(state_ic.get_x_dot(), s_bar_dot_init)
720  );
721 
722  return state_and_sens_ic;
723 
724 }
725 
726 
732 template<class Scalar>
733 inline
734 Thyra::ModelEvaluatorBase::InArgs<Scalar>
736  const ForwardSensitivityStepper<Scalar> &fwdSensStepper,
737  const Thyra::ModelEvaluatorBase::InArgs<Scalar> &state_and_sens_ic
738  )
739 {
740 
741  using Thyra::productVectorBase;
742  typedef Thyra::ModelEvaluatorBase MEB;
743 
744  MEB::InArgs<Scalar>
745  state_ic = fwdSensStepper.getStateModel()->createInArgs();
746 
747  // Copy time, parameters etc.
748  state_ic.setArgs(state_and_sens_ic);
749  state_ic.set_x(
750  productVectorBase(state_and_sens_ic.get_x())->getVectorBlock(0));
751  state_ic.set_x_dot(
752  productVectorBase(state_and_sens_ic.get_x_dot())->getVectorBlock(0));
753 
754  return state_ic;
755 
756 }
757 
758 
759 //
760 // Implementation
761 //
762 
763 
764 // Static members
765 
766 
767 template<class Scalar>
769 = "Force Up-To-Date Jacobian";
770 
771 template<class Scalar>
773 = true;
774 
775 
776 // Constructors, Intializers, Misc.
777 
778 
779 template<class Scalar>
781  :forceUpToDateW_(false),
782  isSingleResidualStepper_(false)
783 {}
784 
785 
786 template<class Scalar>
788  const RCP<const Thyra::ModelEvaluator<Scalar> > &stateModel,
789  const int p_index,
790  const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateBasePoint,
791  const RCP<StepperBase<Scalar> > &stateStepper,
792  const RCP<Thyra::NonlinearSolverBase<Scalar> > &stateTimeStepSolver,
793  const RCP<StepperBase<Scalar> > &sensStepper,
794  const RCP<Thyra::NonlinearSolverBase<Scalar> > &sensTimeStepSolver
795  )
796 
797 {
798  initializeCommon( stateModel, p_index, Teuchos::null, stateBasePoint, stateStepper,
799  stateTimeStepSolver, sensStepper, sensTimeStepSolver );
800 }
801 
802 
803 template<class Scalar>
805  const RCP<const Thyra::ModelEvaluator<Scalar> >& stateModel,
806  const RCP<const Thyra::VectorSpaceBase<Scalar> >& p_space,
807  const Thyra::ModelEvaluatorBase::InArgs<Scalar>& stateBasePoint,
808  const RCP<StepperBase<Scalar> >& stateStepper,
809  const RCP<Thyra::NonlinearSolverBase<Scalar> >& stateTimeStepSolver,
810  const RCP<StepperBase<Scalar> >& sensStepper,
811  const RCP<Thyra::NonlinearSolverBase<Scalar> >& sensTimeStepSolver
812  )
813 {
814  initializeCommon(stateModel, -1, p_space, stateBasePoint, stateStepper,
815  stateTimeStepSolver, sensStepper, sensTimeStepSolver );
816 }
817 
818 
819 template<class Scalar>
821  const RCP<const Thyra::ModelEvaluator<Scalar> > &stateModel,
822  const int p_index,
823  const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateBasePoint,
824  const RCP<StepperBase<Scalar> > &stateStepper,
825  const RCP<Thyra::NonlinearSolverBase<Scalar> > &stateTimeStepSolver,
826  const RCP<IntegratorBase<Scalar> > &stateIntegrator,
827  const Scalar &finalTime,
828  const RCP<StepperBase<Scalar> > &sensStepper,
829  const RCP<Thyra::NonlinearSolverBase<Scalar> > &sensTimeStepSolver
830  )
831 {
832  TEUCHOS_ASSERT(nonnull(stateIntegrator));
833  initializeCommon( stateModel, p_index, Teuchos::null, stateBasePoint, stateStepper,
834  stateTimeStepSolver, sensStepper, sensTimeStepSolver );
835  stateIntegrator_ = stateIntegrator;
836  finalTime_ = finalTime;
837 }
838 
839 
840 template<class Scalar>
842 {
843  return stateModel_.isConst();
844 }
845 
846 
847 template<class Scalar>
848 RCP<const Thyra::ModelEvaluator<Scalar> >
850 {
851  return stateModel_.getConstObj();
852 }
853 
854 
855 template<class Scalar>
856 RCP<StepperBase<Scalar> >
858 {
859  return stateStepper_;
860 }
861 
862 
863 template<class Scalar>
864 RCP<const ForwardSensitivityModelEvaluatorBase<Scalar> >
866 {
867  return sensModel_;
868 }
869 
870 
871 template<class Scalar>
872 RCP<const StateAndForwardSensitivityModelEvaluator<Scalar> >
874 {
875  return stateAndSensModel_;
876 }
877 
878 
879 // Overridden from Teuchos::ParameterListAcceptor
880 
881 
882 template<class Scalar>
884  RCP<Teuchos::ParameterList> const& paramList
885  )
886 {
887  TEUCHOS_TEST_FOR_EXCEPT(is_null(paramList));
888  paramList->validateParameters(*getValidParameters());
889  this->setMyParamList(paramList);
890  forceUpToDateW_ = paramList->get(forceUpToDateW_name_,forceUpToDateW_default_);
891  Teuchos::readVerboseObjectSublist(&*paramList,this);
892 }
893 
894 
895 template<class Scalar>
896 RCP<const Teuchos::ParameterList>
898 {
899  static RCP<const ParameterList> validPL;
900  if (is_null(validPL) ) {
901  RCP<ParameterList> pl = Teuchos::parameterList();
902  pl->set( forceUpToDateW_name_, forceUpToDateW_default_,
903  "If set to true, then the Jacobian matrix W used in the\n"
904  "state timestep equation will be forced to be up to date\n"
905  "with the final value for x for the nonlinear solve. If\n"
906  "you are willing to live with slightly less accurate sensitivities\n"
907  "then set this to false."
908  );
909  Teuchos::setupVerboseObjectSublist(&*pl);
910  validPL = pl;
911  }
912  return validPL;
913 }
914 
915 
916 // Overridden from StepperBase
917 
918 template<class Scalar>
920 {
921  return false;
922 }
923 
924 template<class Scalar>
926  const RCP<const Thyra::ModelEvaluator<Scalar> >& /* model */
927  )
928 {
929  TEUCHOS_TEST_FOR_EXCEPT_MSG( true,
930  "Error, this stepper subclass does not accept a model"
931  " as defined by the StepperBase interface!");
932 }
933 
934 
935 template<class Scalar>
937  const RCP<Thyra::ModelEvaluator<Scalar> >& /* model */
938  )
939 {
940  TEUCHOS_TEST_FOR_EXCEPT_MSG( true,
941  "Error, this stepper subclass does not accept a model"
942  " as defined by the StepperBase interface!");
943 }
944 
945 
946 template<class Scalar>
947 RCP<const Thyra::ModelEvaluator<Scalar> >
949 {
950  return stateAndSensModel_;
951 }
952 
953 
954 template<class Scalar>
955 RCP<Thyra::ModelEvaluator<Scalar> >
957 {
958  return stateAndSensModel_;
959 }
960 
961 
962 template<class Scalar>
964  const Thyra::ModelEvaluatorBase::InArgs<Scalar> &state_and_sens_ic
965  )
966 {
967 
968  typedef Thyra::ModelEvaluatorBase MEB;
969 
970  stateAndSensBasePoint_ = state_and_sens_ic;
971 
972  // Get the product vectors for x_bar = [ x; s_bar ] and x_bar_dot
973 
974  TEUCHOS_TEST_FOR_EXCEPTION(
975  is_null(state_and_sens_ic.get_x()), std::logic_error,
976  "Error, the initial condition for x_bar = [ x; s_bar ] can not be null!" );
977 
978  const RCP<const Thyra::ProductVectorBase<Scalar> >
979  x_bar_init = Thyra::productVectorBase<Scalar>(
980  state_and_sens_ic.get_x()
981  );
982 
983  RCP<const Thyra::ProductVectorBase<Scalar> > x_bar_dot_init;
984  if (state_and_sens_ic.supports(MEB::IN_ARG_x_dot)) {
985  x_bar_dot_init = Thyra::productVectorBase<Scalar>(
986  state_and_sens_ic.get_x_dot()
987  );
988  }
989 
990  // Remove x and x_dot from state_and_sens_ic_in to avoid cloning x and x dot!
991 
992  Thyra::ModelEvaluatorBase::InArgs<Scalar>
993  state_and_sens_ic_no_x = state_and_sens_ic;
994  state_and_sens_ic_no_x.set_x(Teuchos::null);
995  if (state_and_sens_ic_no_x.supports(MEB::IN_ARG_x_dot)) {
996  state_and_sens_ic_no_x.set_x_dot(Teuchos::null);
997  }
998 
999  // Set initial condition for the state
1000 
1001  MEB::InArgs<Scalar> state_ic = stateModel_->createInArgs();
1002  state_ic.setArgs(state_and_sens_ic_no_x,true,true); // Set time, parameters etc.
1003  state_ic.set_x(x_bar_init->getVectorBlock(0)->clone_v());
1004  if (state_ic.supports(MEB::IN_ARG_x_dot)) {
1005  state_ic.set_x_dot(
1006  !is_null(x_bar_dot_init)
1007  ? x_bar_dot_init->getVectorBlock(0)->clone_v()
1008  : Teuchos::null
1009  );
1010  }
1011  stateStepper_->setInitialCondition(state_ic);
1012 
1013  // Set up the integrator if needed
1014  //if (!is_null(stateIntegrator_)) {
1015  // stateIntegrator_->setStepper( stateStepper_, finalTime_ );
1016  // sensModel_->setStateIntegrator( stateIntegrator_, state_ic );
1017  //}
1018 
1019  // Set initial condition for the sensitivities
1020 
1021  MEB::InArgs<Scalar> sens_ic = sensModel_->createInArgs();
1022  sens_ic.setArgs(state_and_sens_ic_no_x,true,true); // Set time etc.
1023  sens_ic.set_x(x_bar_init->getVectorBlock(1)->clone_v());
1024  if (sens_ic.supports(MEB::IN_ARG_x_dot)) {
1025  sens_ic.set_x_dot(
1026  !is_null(x_bar_dot_init)
1027  ? x_bar_dot_init->getVectorBlock(1)->clone_v()
1028  : Teuchos::null
1029  );
1030  }
1031  sensStepper_->setInitialCondition(sens_ic);
1032 
1033 }
1034 
1035 
1036 template<class Scalar>
1037 Thyra::ModelEvaluatorBase::InArgs<Scalar>
1039 {
1040  return stateAndSensBasePoint_;
1041 }
1042 
1043 
1044 template<class Scalar>
1045 Scalar
1047  Scalar dt, StepSizeType stepType
1048  )
1049 {
1050 
1051  RYTHMOS_FUNC_TIME_MONITOR("Rythmos:ForwardSensitivityStepper::takeStep");
1052 
1053  if (!is_null(stateIntegrator_)) {
1054  return takeDecoupledStep(dt,stepType);
1055  }
1056 
1057  return takeSyncedStep(dt,stepType);
1058 
1059 }
1060 
1061 
1062 template<class Scalar>
1063 const StepStatus<Scalar>
1065 {
1066 
1067  const StepStatus<Scalar> sensStepStatus = sensStepper_->getStepStatus();
1068  StepStatus<Scalar> stepStatus;
1069 
1070  stepStatus.message = sensStepStatus.message;
1071  stepStatus.stepStatus = sensStepStatus.stepStatus;
1072  stepStatus.stepLETStatus = sensStepStatus.stepLETStatus;
1073  stepStatus.stepSize = sensStepStatus.stepSize;
1074  stepStatus.order = sensStepStatus.order;
1075  stepStatus.time = sensStepStatus.time;
1076  stepStatus.stepLETValue = sensStepStatus.stepLETValue;
1077  stepStatus.extraParameters = sensStepStatus.extraParameters;
1078 
1079  if (is_null(stateIntegrator_)) {
1080  const StepStatus<Scalar>
1081  stateStepStatus = stateStepper_->getStepStatus();
1082  if (!is_null(stateStepStatus.solution) && !is_null(sensStepStatus.solution))
1083  stepStatus.solution = stateAndSensModel_->create_x_bar_vec(
1084  stateStepStatus.solution, sensStepStatus.solution
1085  );
1086  if (!is_null(stateStepStatus.solutionDot) && !is_null(sensStepStatus.solutionDot))
1087  stepStatus.solutionDot = stateAndSensModel_->create_x_bar_vec(
1088  stateStepStatus.solutionDot, sensStepStatus.solutionDot
1089  );
1090  }
1091 
1092  return stepStatus;
1093 
1094 }
1095 
1096 
1097 // Overridden from InterpolationBufferBase
1098 
1099 
1100 template<class Scalar>
1101 RCP<const Thyra::VectorSpaceBase<Scalar> >
1103 {
1104  return stateAndSensModel_->get_x_space();
1105 }
1106 
1107 
1108 template<class Scalar>
1110  const Array<Scalar>& /* time_vec */,
1111  const Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >& /* x_vec */,
1112  const Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >& /* xdot_vec */
1113  )
1114 {
1115  TEUCHOS_TEST_FOR_EXCEPT("Not implemented addPoints(...) yet but we could if we wanted!");
1116 }
1117 
1118 
1119 template<class Scalar>
1122 {
1123  return sensStepper_->getTimeRange();
1124 }
1125 
1126 
1127 template<class Scalar>
1129  const Array<Scalar>& time_vec,
1130  Array<RCP<const Thyra::VectorBase<Scalar> > >* x_bar_vec,
1131  Array<RCP<const Thyra::VectorBase<Scalar> > >* x_bar_dot_vec,
1132  Array<ScalarMag>* accuracy_vec
1133  ) const
1134 {
1135 
1136  using Teuchos::as;
1137 
1138 #ifdef HAVE_RYTHMOS_DEBUG
1139  TEUCHOS_TEST_FOR_EXCEPT( as<int>(time_vec.size()) == 0 );
1140 #endif
1141 
1142  const int numTimePoints = time_vec.size();
1143 
1144  if (x_bar_vec)
1145  x_bar_vec->clear();
1146 
1147  if (x_bar_dot_vec)
1148  x_bar_dot_vec->clear();
1149 
1150  Array<RCP<const Thyra::VectorBase<Scalar> > >
1151  x_vec, x_dot_vec;
1152 
1153  if (!is_null(stateIntegrator_)) {
1154  stateIntegrator_->getPoints(
1155  time_vec,
1156  x_bar_vec ? &x_vec: 0,
1157  x_bar_dot_vec ? &x_dot_vec: 0,
1158  0 // Ignoring accuracy from state for now!
1159  );
1160  }
1161  else {
1162  stateStepper_->getPoints(
1163  time_vec,
1164  x_bar_vec ? &x_vec: 0,
1165  x_bar_dot_vec ? &x_dot_vec: 0,
1166  0 // Ignoring accuracy from state for now!
1167  );
1168  }
1169 
1170  Array<RCP<const Thyra::VectorBase<Scalar> > >
1171  s_bar_vec, s_bar_dot_vec;
1172 
1173  sensStepper_->getPoints(
1174  time_vec,
1175  x_bar_vec ? &s_bar_vec: 0,
1176  x_bar_dot_vec ? &s_bar_dot_vec: 0,
1177  accuracy_vec
1178  );
1179 
1180  if ( x_bar_vec ) {
1181  for ( int i = 0; i < numTimePoints; ++i ) {
1182  x_bar_vec->push_back(
1183  stateAndSensModel_->create_x_bar_vec(x_vec[i],s_bar_vec[i])
1184  );
1185  }
1186  }
1187 
1188  if ( x_bar_dot_vec ) {
1189  for ( int i = 0; i < numTimePoints; ++i ) {
1190  x_bar_dot_vec->push_back(
1191  stateAndSensModel_->create_x_bar_vec(x_dot_vec[i],s_bar_dot_vec[i])
1192  );
1193  }
1194  }
1195 
1196 }
1197 
1198 
1199 template<class Scalar>
1201  Array<Scalar>* time_vec
1202  ) const
1203 {
1204  TEUCHOS_ASSERT( time_vec != NULL );
1205  time_vec->clear();
1206  if (is_null(stateIntegrator_) && is_null(stateStepper_)) {
1207  return;
1208  }
1209  if (!is_null(stateIntegrator_)) {
1210  stateIntegrator_->getNodes(time_vec);
1211  }
1212  else {
1213  stateStepper_->getNodes(time_vec);
1214  }
1215 }
1216 
1217 
1218 template<class Scalar>
1220  Array<Scalar>& /* time_vec */
1221  )
1222 {
1223  TEUCHOS_TEST_FOR_EXCEPT("Not implemented yet but we can!");
1224 }
1225 
1226 
1227 template<class Scalar>
1229 {
1230  return sensStepper_->getOrder();
1231  // Note: This assumes that stateStepper will have the same order!
1232 }
1233 
1234 
1235 // private
1236 
1237 
1238 template<class Scalar>
1240  const RCP<const Thyra::ModelEvaluator<Scalar> >& stateModel,
1241  const int p_index,
1242  const RCP<const Thyra::VectorSpaceBase<Scalar> > &p_space,
1243  const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateBasePoint,
1244  const RCP<StepperBase<Scalar> > &stateStepper,
1245  const RCP<Thyra::NonlinearSolverBase<Scalar> > &stateTimeStepSolver,
1246  const RCP<StepperBase<Scalar> > &sensStepper,
1247  const RCP<Thyra::NonlinearSolverBase<Scalar> > &sensTimeStepSolver
1248  )
1249 {
1250 
1251  using Teuchos::rcp_implicit_cast;
1252  using Teuchos::rcp_dynamic_cast;
1253 
1254  typedef Thyra::ModelEvaluatorBase MEB;
1255 
1256  //
1257  // Validate input
1258  //
1259 
1260  TEUCHOS_ASSERT( p_index >= 0 || nonnull(p_space) );
1261  if (nonnull(p_space)) {
1262  TEUCHOS_ASSERT_EQUALITY(p_index, -1);
1263  }
1264  if (p_index >= 0) {
1265  TEUCHOS_ASSERT(is_null(p_space));
1266  }
1267  TEUCHOS_TEST_FOR_EXCEPT( is_null(stateModel) );
1268  TEUCHOS_TEST_FOR_EXCEPT( is_null(stateStepper) );
1269  if (stateStepper->isImplicit()) {
1270  TEUCHOS_TEST_FOR_EXCEPT( is_null(stateTimeStepSolver) ); // allow to be null for explicit methods
1271  }
1272 
1273  //
1274  // Create the sensModel which will do some more validation
1275  //
1276 
1277  RCP<ForwardSensitivityModelEvaluatorBase<Scalar> > sensModel;
1278  MEB::InArgs<Scalar> stateModelInArgs = stateModel->createInArgs();
1279  if (stateModelInArgs.supports(MEB::IN_ARG_x_dot)) {
1280  // Implicit DE formulation
1281  sensModel = Teuchos::rcp(new ForwardSensitivityImplicitModelEvaluator<Scalar>);
1282  }
1283  else {
1284  // Explicit DE formulation
1285  sensModel = Teuchos::rcp(new ForwardSensitivityExplicitModelEvaluator<Scalar>);
1286  }
1287 
1288  if (p_index >= 0) {
1289  sensModel->initializeStructure(stateModel, p_index);
1290  }
1291  else {
1292  sensModel->initializeStructureInitCondOnly(stateModel, p_space);
1293  }
1294 
1295  //
1296  // Get the input objects
1297  //
1298 
1299  stateModel_.initialize(stateModel);
1300 
1301  stateBasePoint_ = stateBasePoint;
1302 
1303  stateStepper_ = stateStepper;
1304 
1305  stateTimeStepSolver_ = stateTimeStepSolver;
1306 
1307  sensModel_ = sensModel;
1308 
1309  stateAndSensModel_ = Teuchos::rcp(new StateAndForwardSensitivityModelEvaluator<Scalar>);
1310  stateAndSensModel_->initializeStructure(sensModel_);
1311 
1312  if (!is_null(sensStepper)) {
1313  sensStepper_ = sensStepper;
1314  }
1315  else {
1316  sensStepper_ = stateStepper_->cloneStepperAlgorithm();
1317  TEUCHOS_TEST_FOR_EXCEPTION(
1318  is_null(sensStepper_), std::logic_error,
1319  "Error, if the client does not pass in a stepper for the senitivity\n"
1320  "equations then the stateStepper object must support cloning to create\n"
1321  "the sensitivity stepper!"
1322  );
1323  }
1324 
1325  if (!is_null(sensTimeStepSolver)) {
1326  sensTimeStepSolver_ = sensTimeStepSolver;
1327  }
1328  else {
1329  RCP<Thyra::LinearNonlinearSolver<Scalar> >
1330  linearNonlinearSolver(new Thyra::LinearNonlinearSolver<Scalar>);
1331  // ToDo: Set tolerance on the nonlinear solver???
1332  sensTimeStepSolver_ = linearNonlinearSolver;
1333  }
1334 
1335  //
1336  // Setup the steppers
1337  //
1338 
1339  isSingleResidualStepper_ = true; // ToDo: Add dynamic cast on
1340  // stateTimeStepSolver to check this!
1341 
1342  setStepperModel(Teuchos::inOutArg(*stateStepper_),stateModel_);
1343  if (stateStepper_->isImplicit()) {
1344  rcp_dynamic_cast<SolverAcceptingStepperBase<Scalar> >(
1345  stateStepper_,true)->setSolver(stateTimeStepSolver_);
1346  }
1347  sensStepper_->setModel(sensModel_);
1348  if (sensStepper_->isImplicit()) {
1349  rcp_dynamic_cast<SolverAcceptingStepperBase<Scalar> >(
1350  sensStepper_,true)->setSolver(sensTimeStepSolver_);
1351  }
1352 
1353  stateBasePoint_t_ = stateModel_->createInArgs();
1354 
1355  // 2007/05/18: rabartl: ToDo: Move the above initialization code to give
1356  // setInitializeCondition(...) a chance to set the initial condition.
1357 
1358 }
1359 
1360 
1361 template<class Scalar>
1362 Scalar ForwardSensitivityStepper<Scalar>::takeSyncedStep(
1363  Scalar dt, StepSizeType stepType
1364  )
1365 {
1366 
1367  RYTHMOS_FUNC_TIME_MONITOR_DIFF("Rythmos:ForwardSensitivityStepper::takeStep: synced",
1368  TopLevel);
1369 
1370  using Teuchos::as;
1371  typedef Teuchos::ScalarTraits<Scalar> ST;
1372  typedef Teuchos::VerboseObjectTempState<InterpolationBufferBase<Scalar> > VOTSIBB;
1373 
1374  RCP<Teuchos::FancyOStream> out = this->getOStream();
1375  Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
1376  const bool lowTrace =
1377  ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) );
1378  const bool mediumTrace =
1379  ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_MEDIUM) );
1380  Teuchos::OSTab tab(out);
1381 
1382  if (lowTrace) {
1383  *out
1384  << "\nEntering " << TypeNameTraits<ForwardSensitivityStepper<Scalar> >::name()
1385  << "::takeSyncedStep("<<dt<<","<<toString(stepType)<<") ...\n";
1386  }
1387 
1388  //
1389  // A) Compute the state timestep
1390  //
1391 
1392  if (lowTrace) {
1393  *out
1394  << "\nTaking state step using stepper : "
1395  << stateStepper_->description() << "\n";
1396  }
1397 
1398  Scalar state_dt = -1.0;
1399  {
1400  RYTHMOS_FUNC_TIME_MONITOR("Rythmos:ForwardSensitivityStepper::takeStep: stateStep");
1401  VOTSIBB stateStepper_outputTempState(stateStepper_,out,verbLevel);
1402  state_dt = stateStepper_->takeStep(dt,stepType);
1403  }
1404 
1405  if (state_dt < Scalar(-ST::one())) {
1406  if (lowTrace)
1407  *out << "\nThe state stepper has failed so return a failed timestep!\n";
1408  return state_dt;
1409  }
1410 
1411  {
1412  RYTHMOS_FUNC_TIME_MONITOR("Rythmos:ForwardSensitivityStepper::takeStep: updateSensModel");
1413  // Set up the sensitivity model for this timestep
1414  sensModel_->initializePointState(Teuchos::inOutArg(*stateStepper_),forceUpToDateW_);
1415  }
1416 
1417  //
1418  // C) Compute the sensitivity timestep for the exact same timestep as was
1419  // used for the state solve.
1420  //
1421 
1422  if (lowTrace) {
1423  *out
1424  << "\nTaking sensitivity step using stepper : "
1425  << sensStepper_->description() << "\n";
1426  }
1427 
1428  Scalar sens_dt = -1.0;
1429  {
1430  RYTHMOS_FUNC_TIME_MONITOR("Rythmos:ForwardSensitivityStepper::takeStep: sensStep");
1431  // Copy the step control data to make sure that the sensStepper takes the
1432  // same type of step that the statStepper took. This is needed to ensure
1433  // that the W matrix is the same for one.
1434  sensStepper_->setStepControlData(*stateStepper_);
1435  VOTSIBB sensStepper_outputTempState(sensStepper_,out,verbLevel);
1436  sens_dt = sensStepper_->takeStep(state_dt,STEP_TYPE_FIXED);
1437  }
1438 
1439  if (mediumTrace) {
1440  const StepStatus<Scalar> sensStepStatus = sensStepper_->getStepStatus();
1441  *out << "\nSensitivity step status:\n" << sensStepStatus;
1442  }
1443 
1444  TEUCHOS_TEST_FOR_EXCEPTION(
1445  sens_dt != state_dt, std::logic_error,
1446  "Error, the sensitivity step failed for some reason. We should\n"
1447  "just return a negative step size and reject the step but currently\n"
1448  "there is no way to roll back the state timestep it for back to\n"
1449  "the status before this function was called!"
1450  );
1451 
1452  // 2007/05/18: rabartl: ToDo: If stepType == STEP_TYPE_VARIABLE and the state
1453  // timestep sucessed but the sensitivity timestep failed, then we need to
1454  // really throw an excpetion because there is nothing that we can really do
1455  // here!
1456 
1457  // 2007/05/18: rabartl: ToDo: Replace the above std::logic_error type with
1458  // a Rythmos::CatastrophicFailure or just use Thyra::CatastrophicFailure!
1459 
1460  if (lowTrace) {
1461  *out
1462  << "\nLeaving " << TypeNameTraits<ForwardSensitivityStepper<Scalar> >::name()
1463  << "::takeSyncedStep("<<dt<<","<<toString(stepType)<<") ...\n";
1464  }
1465 
1466  return state_dt;
1467 
1468 }
1469 
1470 
1471 template<class Scalar>
1472 Scalar ForwardSensitivityStepper<Scalar>::takeDecoupledStep(
1473  Scalar dt, StepSizeType stepType
1474  )
1475 {
1476 
1477  RYTHMOS_FUNC_TIME_MONITOR("Rythmos:ForwardSensitivityStepper::takeStep: decoupled");
1478 
1479  using Teuchos::as;
1480  typedef Teuchos::VerboseObjectTempState<InterpolationBufferBase<Scalar> > VOTSIBB;
1481 
1482  RCP<Teuchos::FancyOStream> out = this->getOStream();
1483  Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
1484  const bool lowTrace =
1485  ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) );
1486  const bool mediumTrace =
1487  ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_MEDIUM) );
1488  Teuchos::OSTab tab(out);
1489 
1490  if (lowTrace) {
1491  *out
1492  << "\nEntering " << TypeNameTraits<ForwardSensitivityStepper<Scalar> >::name()
1493  << "::takeDecoupledStep("<<dt<<","<<toString(stepType)<<") ...\n";
1494  }
1495 
1496  //
1497  // A) Take the sens timestep
1498  //
1499 
1500  if (lowTrace) {
1501  *out
1502  << "\nTaking sensitivity step using stepper : "
1503  << sensStepper_->description() << "\n";
1504  }
1505 
1506  Scalar sens_dt = -1.0;
1507  VOTSIBB sensStepper_outputTempState(sensStepper_,out,verbLevel);
1508  sens_dt = sensStepper_->takeStep(dt,stepType);
1509 
1510  if (mediumTrace) {
1511  const StepStatus<Scalar> sensStepStatus = sensStepper_->getStepStatus();
1512  *out << "\nSensitivity step status:\n" << sensStepStatus;
1513  }
1514 
1515  //
1516  // B) Wipe out all state interp buffer info before this sens timestep
1517  //
1518 
1519  //TEUCHOS_TEST_FOR_EXCEPT(true);
1520 
1521  if (lowTrace) {
1522  *out
1523  << "\nLeaving " << TypeNameTraits<ForwardSensitivityStepper<Scalar> >::name()
1524  << "::takeDecoupledStep("<<dt<<","<<toString(stepType)<<") ...\n";
1525  }
1526 
1527  return sens_dt;
1528 
1529 }
1530 
1531 
1532 } // namespace Rythmos
1533 
1534 
1535 #endif //RYTHMOS_FORWARD_SENSITIVITY_STEPPER_HPP
void initialize(const RCP< const Thyra::ModelEvaluator< Scalar > > &stateModel, const int p_index, const Thyra::ModelEvaluatorBase::InArgs< Scalar > &stateBasePoint, const RCP< StepperBase< Scalar > > &stateStepper, const RCP< Thyra::NonlinearSolverBase< Scalar > > &stateTimeStepSolver, const RCP< StepperBase< Scalar > > &sensStepper=Teuchos::null, const RCP< Thyra::NonlinearSolverBase< Scalar > > &sensTimeStepSolver=Teuchos::null)
Deprecated.
RCP< Thyra::ModelEvaluator< Scalar > > getNonconstModel()
RCP< const Thyra::ModelEvaluator< Scalar > > getModel() const
Returns getStateAndFwdSensModel().
const StepStatus< Scalar > getStepStatus() const
Base class for defining stepper functionality.
Abstract interface for time integrators.
RCP< StepperBase< Scalar > > getNonconstStateStepper()
Return the state stepper that was passed into an initialize function.
RCP< const Teuchos::ParameterList > extraParameters
bool stateModelIsConst() const
Return if the state model is const-only or not.
Foward sensitivity stepper concrete subclass.
void setInitialCondition(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &state_and_sens_ic)
Sets the full initial condition for x_bar = [ x; s_bar] and x_bar_dot = [ x_dot; s_bar_dot ] as well...
void initializeSyncedSteppers(const RCP< const Thyra::ModelEvaluator< Scalar > > &stateModel, const int p_index, const Thyra::ModelEvaluatorBase::InArgs< Scalar > &stateBasePoint, const RCP< StepperBase< Scalar > > &stateStepper, const RCP< Thyra::NonlinearSolverBase< Scalar > > &stateTimeStepSolver, const RCP< StepperBase< Scalar > > &sensStepper=Teuchos::null, const RCP< Thyra::NonlinearSolverBase< Scalar > > &sensTimeStepSolver=Teuchos::null)
Intialize for synced state and sens steppers.
RCP< ForwardSensitivityStepper< Scalar > > forwardSensitivityStepper()
Nonmember constructor.
RCP< const ForwardSensitivityModelEvaluatorBase< Scalar > > getFwdSensModel() const
Return the forward sensitivity model evaluator object that got created internally when the initialize...
RCP< ForwardSensitivityStepper< Scalar > > forwardSensitivityStepper(const RCP< const Thyra::ModelEvaluator< Scalar > > &stateModel, const int p_index, const Thyra::ModelEvaluatorBase::InArgs< Scalar > &stateBasePoint, const RCP< StepperBase< Scalar > > &stateStepper, const RCP< Thyra::NonlinearSolverBase< Scalar > > &stateTimeStepSolver, const RCP< StepperBase< Scalar > > &sensStepper=Teuchos::null, const RCP< Thyra::NonlinearSolverBase< Scalar > > &sensTimeStepSolver=Teuchos::null)
Nonmember constructor.
void initializeSyncedSteppersInitCondOnly(const RCP< const Thyra::ModelEvaluator< Scalar > > &stateModel, const RCP< const Thyra::VectorSpaceBase< Scalar > > &p_space, const Thyra::ModelEvaluatorBase::InArgs< Scalar > &stateBasePoint, const RCP< StepperBase< Scalar > > &stateStepper, const RCP< Thyra::NonlinearSolverBase< Scalar > > &stateTimeStepSolver, const RCP< StepperBase< Scalar > > &sensStepper=Teuchos::null, const RCP< Thyra::NonlinearSolverBase< Scalar > > &sensTimeStepSolver=Teuchos::null)
Intialize for synced state and sens steppers for an initial-condition only parametrized sensitivity p...
RCP< const Thyra::VectorBase< Scalar > > solutionDot
Thyra::ModelEvaluatorBase::InArgs< Scalar > createStateAndSensInitialCondition(const ForwardSensitivityStepper< Scalar > &fwdSensStepper, const Thyra::ModelEvaluatorBase::InArgs< Scalar > &state_ic, const RCP< const Thyra::MultiVectorBase< Scalar > > S_init=Teuchos::null, const RCP< const Thyra::MultiVectorBase< Scalar > > S_dot_init=Teuchos::null)
Set up default initial conditions for the state and sensitivity stepper with default zero initial con...
Thyra::ModelEvaluatorBase::InArgs< Scalar > extractStateInitialCondition(const ForwardSensitivityStepper< Scalar > &fwdSensStepper, const Thyra::ModelEvaluatorBase::InArgs< Scalar > &state_and_sens_ic)
Extract out the initial condition for just the state model given the initial condition for the state ...
void setModel(const RCP< const Thyra::ModelEvaluator< Scalar > > &model)
Throws exception.
RCP< const Teuchos::ParameterList > getValidParameters() const
RCP< const Thyra::ModelEvaluator< Scalar > > getStateModel() const
Return the state model that was passed into an initialize function.
ForwardSensitivityStepper()
Constructs to uninitialized.
Scalar takeStep(Scalar dt, StepSizeType stepType)
Teuchos::ScalarTraits< Scalar >::magnitudeType ScalarMag
RCP< const Thyra::VectorBase< Scalar > > solution
void setParameterList(RCP< Teuchos::ParameterList > const &paramList)
void getNodes(Array< Scalar > *time_vec) const
RCP< const Thyra::VectorSpaceBase< Scalar > > get_x_space() const
Returns the space for x_bar and x_bar_dot.
void setNonconstModel(const RCP< Thyra::ModelEvaluator< Scalar > > &model)
Throws exception.
RCP< const StateAndForwardSensitivityModelEvaluator< Scalar > > getStateAndFwdSensModel() const
Return the state and forward sensitivity model evaluator object that got created internally when the ...
void addPoints(const Array< Scalar > &time_vec, const Array< Teuchos::RCP< const Thyra::VectorBase< Scalar > > > &x_vec, const Array< Teuchos::RCP< const Thyra::VectorBase< Scalar > > > &xdot_vec)
int getParameterIndex(const ForwardSensitivityStepper< Scalar > &fwdSensStepper)
Return the index of the parameter subvector in the underlying state model.
void getPoints(const Array< Scalar > &time_vec, Array< RCP< const Thyra::VectorBase< Scalar > > > *x_vec, Array< RCP< const Thyra::VectorBase< Scalar > > > *xdot_vec, Array< ScalarMag > *accuracy_vec) const
Thyra::ModelEvaluatorBase::InArgs< Scalar > getInitialCondition() const
void initializeDecoupledSteppers(const RCP< const Thyra::ModelEvaluator< Scalar > > &stateModel, const int p_index, const Thyra::ModelEvaluatorBase::InArgs< Scalar > &stateBasePoint, const RCP< StepperBase< Scalar > > &stateStepper, const RCP< Thyra::NonlinearSolverBase< Scalar > > &stateTimeStepSolver, const RCP< IntegratorBase< Scalar > > &stateIntegrator, const Scalar &finalTime, const RCP< StepperBase< Scalar > > &sensStepper=Teuchos::null, const RCP< Thyra::NonlinearSolverBase< Scalar > > &sensTimeStepSolver=Teuchos::null)
Intialize for decoupled state and sens steppers.