Rythmos - Transient Integration for Differential Equations Version of the Day
Rythmos_ForwardSensitivityStepper.hpp
00001 //@HEADER
00002 // ***********************************************************************
00003 //
00004 //                           Rythmos Package
00005 //                 Copyright (2006) Sandia Corporation
00006 //
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Todd S. Coffey (tscoffe@sandia.gov)
00025 //
00026 // ***********************************************************************
00027 //@HEADER
00028 
00029 #ifndef RYTHMOS_FORWARD_SENSITIVITY_STEPPER_HPP
00030 #define RYTHMOS_FORWARD_SENSITIVITY_STEPPER_HPP
00031 
00032 
00033 #include "Rythmos_StepperBase.hpp"
00034 #include "Rythmos_StepperHelpers.hpp"
00035 #include "Rythmos_ForwardSensitivityModelEvaluatorBase.hpp"
00036 #include "Rythmos_ForwardSensitivityImplicitModelEvaluator.hpp"
00037 #include "Rythmos_ForwardSensitivityExplicitModelEvaluator.hpp"
00038 #include "Rythmos_StateAndForwardSensitivityModelEvaluator.hpp"
00039 #include "Rythmos_SolverAcceptingStepperBase.hpp"
00040 #include "Rythmos_IntegratorBase.hpp"
00041 #include "Rythmos_SingleResidualModelEvaluatorBase.hpp"
00042 #include "Thyra_ModelEvaluatorHelpers.hpp"
00043 #include "Thyra_LinearNonlinearSolver.hpp"
00044 #include "Thyra_ProductVectorBase.hpp"
00045 #include "Thyra_AssertOp.hpp"
00046 #include "Teuchos_ParameterListAcceptorDefaultBase.hpp"
00047 #include "Teuchos_VerboseObjectParameterListHelpers.hpp"
00048 #include "Teuchos_ConstNonconstObjectContainer.hpp"
00049 #include "Teuchos_Assert.hpp"
00050 #include "Teuchos_as.hpp"
00051 
00052 
00053 namespace Rythmos {
00054 
00055 
00208 template<class Scalar>
00209 class ForwardSensitivityStepper
00210   : virtual public StepperBase<Scalar>,
00211     virtual public Teuchos::ParameterListAcceptorDefaultBase
00212 {
00213 public:
00214 
00216   typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType ScalarMag;
00217 
00220 
00222   ForwardSensitivityStepper();
00223 
00277   void initializeSyncedSteppers(
00278     const RCP<const Thyra::ModelEvaluator<Scalar> > &stateModel,
00279     const int p_index,
00280     const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateBasePoint,
00281     const RCP<StepperBase<Scalar> > &stateStepper,
00282     const RCP<Thyra::NonlinearSolverBase<Scalar> > &stateTimeStepSolver,
00283     const RCP<StepperBase<Scalar> > &sensStepper = Teuchos::null,
00284     const RCP<Thyra::NonlinearSolverBase<Scalar> > &sensTimeStepSolver = Teuchos::null
00285     );
00286 
00311   void initializeSyncedSteppersInitCondOnly(
00312     const RCP<const Thyra::ModelEvaluator<Scalar> >& stateModel,
00313     const RCP<const Thyra::VectorSpaceBase<Scalar> >& p_space,
00314     const Thyra::ModelEvaluatorBase::InArgs<Scalar>& stateBasePoint,
00315     const RCP<StepperBase<Scalar> >& stateStepper,
00316     const RCP<Thyra::NonlinearSolverBase<Scalar> >& stateTimeStepSolver,
00317     const RCP<StepperBase<Scalar> >& sensStepper = Teuchos::null,
00318     const RCP<Thyra::NonlinearSolverBase<Scalar> >& sensTimeStepSolver = Teuchos::null
00319     );
00320 
00356   void initializeDecoupledSteppers(
00357     const RCP<const Thyra::ModelEvaluator<Scalar> > &stateModel,
00358     const int p_index,
00359     const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateBasePoint,
00360     const RCP<StepperBase<Scalar> > &stateStepper,
00361     const RCP<Thyra::NonlinearSolverBase<Scalar> > &stateTimeStepSolver,
00362     const RCP<IntegratorBase<Scalar> > &stateIntegrator,
00363     const Scalar &finalTime,
00364     const RCP<StepperBase<Scalar> > &sensStepper = Teuchos::null,
00365     const RCP<Thyra::NonlinearSolverBase<Scalar> > &sensTimeStepSolver = Teuchos::null
00366     );
00367 
00369   bool stateModelIsConst() const;
00370 
00374   RCP<const Thyra::ModelEvaluator<Scalar> >
00375   getStateModel() const;
00376 
00380   RCP<StepperBase<Scalar> >
00381   getNonconstStateStepper();
00382 
00386   RCP<const ForwardSensitivityModelEvaluatorBase<Scalar> >
00387   getFwdSensModel() const;
00388 
00395   RCP<const StateAndForwardSensitivityModelEvaluator<Scalar> >
00396   getStateAndFwdSensModel() const;
00397 
00399 
00402 
00404   void setParameterList(RCP<Teuchos::ParameterList> const& paramList);
00406   RCP<const Teuchos::ParameterList> getValidParameters() const;
00407 
00409 
00412 
00414   bool acceptsModel() const;
00415 
00417   void setModel(
00418     const RCP<const Thyra::ModelEvaluator<Scalar> >& model
00419     );
00420 
00422   void setNonconstModel(
00423     const RCP<Thyra::ModelEvaluator<Scalar> >& model
00424     );
00425 
00432   RCP<const Thyra::ModelEvaluator<Scalar> > getModel() const;
00433 
00435   RCP<Thyra::ModelEvaluator<Scalar> > getNonconstModel();
00436 
00437   // RAB: ToDo: 2007/05/15: I need to talk with Todd about potentially
00438   // removing the setModel() and getModel() functions from the StepperBase
00439   // interface.  In the case of this forward sensitivity solver, I am not sure
00440   // that it makes a lot of sense to define a model.  This surely will not be
00441   // the model that a generic client would expect.  The assumption I am sure
00442   // would be that this model has the same space for x as the interpolation
00443   // buffer but that is not true in this case.
00444 
00458   void setInitialCondition(
00459     const Thyra::ModelEvaluatorBase::InArgs<Scalar> &state_and_sens_ic
00460     );
00461 
00463   Thyra::ModelEvaluatorBase::InArgs<Scalar> getInitialCondition() const;
00464 
00466   Scalar takeStep( Scalar dt, StepSizeType stepType );
00467 
00469   const StepStatus<Scalar> getStepStatus() const;
00470 
00472 
00475 
00482   RCP<const Thyra::VectorSpaceBase<Scalar> >
00483   get_x_space() const;
00484 
00486   void addPoints(
00487     const Array<Scalar>& time_vec,
00488     const Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >& x_vec,
00489     const Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >& xdot_vec
00490     );
00491 
00493   TimeRange<Scalar> getTimeRange() const;
00494 
00496   void getPoints(
00497     const Array<Scalar>& time_vec,
00498     Array<RCP<const Thyra::VectorBase<Scalar> > >* x_vec,
00499     Array<RCP<const Thyra::VectorBase<Scalar> > >* xdot_vec,
00500     Array<ScalarMag>* accuracy_vec
00501     ) const;
00502 
00504   void getNodes(Array<Scalar>* time_vec) const;
00505 
00507   void removeNodes(Array<Scalar>& time_vec);
00508 
00510   int getOrder() const;
00511 
00513 
00516 
00518   void initialize(
00519     const RCP<const Thyra::ModelEvaluator<Scalar> > &stateModel,
00520     const int p_index,
00521     const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateBasePoint,
00522     const RCP<StepperBase<Scalar> > &stateStepper,
00523     const RCP<Thyra::NonlinearSolverBase<Scalar> > &stateTimeStepSolver,
00524     const RCP<StepperBase<Scalar> > &sensStepper = Teuchos::null,
00525     const RCP<Thyra::NonlinearSolverBase<Scalar> > &sensTimeStepSolver = Teuchos::null
00526     )
00527     {
00528       initializeSyncedSteppers(
00529         stateModel, p_index, stateBasePoint, stateStepper, stateTimeStepSolver,
00530         sensStepper, sensTimeStepSolver
00531         );
00532     }
00533 
00535 
00536 private:
00537   // ///////////////////
00538   // Private types
00539 
00540   typedef Teuchos::ConstNonconstObjectContainer<Thyra::ModelEvaluator<Scalar> > CNCME;
00541 
00542   // /////////////////////////
00543   // Private data members
00544 
00545   bool forceUpToDateW_;
00546   CNCME stateModel_;
00547   Thyra::ModelEvaluatorBase::InArgs<Scalar> stateBasePoint_;
00548   RCP<StepperBase<Scalar> > stateStepper_;
00549   RCP<Thyra::NonlinearSolverBase<Scalar> > stateTimeStepSolver_;
00550   RCP<IntegratorBase<Scalar> > stateIntegrator_;
00551   Scalar finalTime_;
00552   Thyra::ModelEvaluatorBase::InArgs<Scalar> stateAndSensBasePoint_;
00553   RCP<StepperBase<Scalar> > sensStepper_;
00554   RCP<Thyra::NonlinearSolverBase<Scalar> > sensTimeStepSolver_;
00555 
00556   bool isSingleResidualStepper_;
00557   RCP<ForwardSensitivityModelEvaluatorBase<Scalar> > sensModel_;
00558   RCP<StateAndForwardSensitivityModelEvaluator<Scalar> > stateAndSensModel_;
00559   Thyra::ModelEvaluatorBase::InArgs<Scalar> stateBasePoint_t_;
00560 
00561   static const std::string forceUpToDateW_name_;
00562   static const bool forceUpToDateW_default_;
00563 
00564   // /////////////////////////
00565   // Private member functions
00566 
00567 
00568   // Common initialization helper
00569   //
00570   // Preconditions:
00571   // (*) p_index >=0 or nonnull(p_space) == true
00572   //
00573   void initializeCommon(
00574     const RCP<const Thyra::ModelEvaluator<Scalar> > &stateModel,
00575     const int p_index,
00576     const RCP<const Thyra::VectorSpaceBase<Scalar> > &p_space,
00577     const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateBasePoint,
00578     const RCP<StepperBase<Scalar> > &stateStepper,
00579     const RCP<Thyra::NonlinearSolverBase<Scalar> > &stateTimeStepSolver,
00580     const RCP<StepperBase<Scalar> > &sensStepper,
00581     const RCP<Thyra::NonlinearSolverBase<Scalar> > &sensTimeStepSolver
00582     );
00583 
00584   Scalar takeSyncedStep( Scalar dt, StepSizeType stepType );
00585 
00586   Scalar takeDecoupledStep( Scalar dt, StepSizeType stepType );
00587 
00588 };
00589 
00590 
00591 // 2009/09/05: rabartl: ToDo: To fix the const and non-const handling of the
00592 // stateModel in this class is going to be a lot of work but here is what
00593 // needs to be done:
00594 //
00595 // (*) Duplicate each function that sets the stateModel, one for const and one
00596 // for non-const.
00597 //
00598 // (*) Create a single a private version for each of these functions that
00599 // accepts a Teuchos::ConstNonconstObjectContainer<> object and will implement
00600 // the guts of the set up same as the existing functions.
00601 //
00602 // (*) Get all of the concrete StepperBase subclasses to implement the
00603 // setModel(const RCP<const ME>&) and modelIsConst() functions and get them to
00604 // use the Teuchos::ConstNonconstObjectContainer<> class as described above.
00605 // This should be pretty easy as the only function that needs to be addressed
00606 // in most cases is just the setModel(...) function.
00607 //
00608 
00609 
00614 template<class Scalar>
00615 inline
00616 RCP<ForwardSensitivityStepper<Scalar> >
00617 forwardSensitivityStepper()
00618 {
00619   return Teuchos::rcp(new ForwardSensitivityStepper<Scalar>());
00620 }
00621 
00622 
00627 template<class Scalar>
00628 inline
00629 RCP<ForwardSensitivityStepper<Scalar> >
00630 forwardSensitivityStepper(
00631   const RCP<const Thyra::ModelEvaluator<Scalar> > &stateModel,
00632   const int p_index,
00633   const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateBasePoint,
00634   const RCP<StepperBase<Scalar> > &stateStepper,
00635   const RCP<Thyra::NonlinearSolverBase<Scalar> > &stateTimeStepSolver,
00636   const RCP<StepperBase<Scalar> > &sensStepper = Teuchos::null,
00637   const RCP<Thyra::NonlinearSolverBase<Scalar> > &sensTimeStepSolver = Teuchos::null
00638   )
00639 {
00640   RCP<ForwardSensitivityStepper<Scalar> >
00641     fwdSensStepper = Teuchos::rcp(new ForwardSensitivityStepper<Scalar>());
00642   fwdSensStepper->initializeSyncedSteppers(
00643     stateModel, p_index, stateBasePoint, stateStepper, stateTimeStepSolver );
00644   return fwdSensStepper;
00645 }
00646 
00647 
00653 template<class Scalar>
00654 int getParameterIndex(
00655   const ForwardSensitivityStepper<Scalar> &fwdSensStepper
00656   )
00657 {
00658   return fwdSensStepper.getFwdSensModel()->get_p_index();
00659 }
00660 
00661 
00668 template<class Scalar>
00669 inline
00670 Thyra::ModelEvaluatorBase::InArgs<Scalar>
00671 createStateAndSensInitialCondition(
00672   const ForwardSensitivityStepper<Scalar> &fwdSensStepper,
00673   const Thyra::ModelEvaluatorBase::InArgs<Scalar> &state_ic,
00674   const RCP<const Thyra::MultiVectorBase<Scalar> > S_init = Teuchos::null,
00675   const RCP<const Thyra::MultiVectorBase<Scalar> > S_dot_init = Teuchos::null
00676   )
00677 {
00678 
00679   using Teuchos::outArg;
00680   using Thyra::assign;
00681   typedef Thyra::ModelEvaluatorBase MEB;
00682 
00683   RCP<const Thyra::VectorBase<Scalar> > s_bar_init;
00684   if (nonnull(S_init)) {
00685     s_bar_init = create_s_bar_given_S(*fwdSensStepper.getFwdSensModel(), S_init);
00686   }
00687   else {
00688     RCP<Thyra::VectorBase<Scalar> > s_bar_init_loc =
00689       createMember(fwdSensStepper.getFwdSensModel()->get_x_space());
00690     assign( outArg(*s_bar_init_loc), 0.0 );
00691     s_bar_init = s_bar_init_loc;
00692   }
00693 
00694   RCP<const Thyra::VectorBase<Scalar> > s_bar_dot_init;
00695   if (nonnull(S_dot_init)) {
00696     s_bar_dot_init = create_s_bar_given_S(*fwdSensStepper.getFwdSensModel(), S_dot_init);
00697   }
00698   else {
00699     RCP<Thyra::VectorBase<Scalar> > s_bar_dot_init_loc =
00700       createMember(fwdSensStepper.getFwdSensModel()->get_x_space());
00701     assign( outArg(*s_bar_dot_init_loc), 0.0 );
00702     s_bar_dot_init = s_bar_dot_init_loc;
00703   }
00704 
00705   RCP<const Rythmos::StateAndForwardSensitivityModelEvaluator<Scalar> >
00706     stateAndSensModel = fwdSensStepper.getStateAndFwdSensModel();
00707 
00708   MEB::InArgs<Scalar>
00709     state_and_sens_ic = fwdSensStepper.getModel()->createInArgs();
00710 
00711   // Copy time, parameters etc.
00712   state_and_sens_ic.setArgs(state_ic);
00713   // Set initial condition for x_bar = [ x; s_bar ]
00714   state_and_sens_ic.set_x(
00715     stateAndSensModel->create_x_bar_vec(state_ic.get_x(), s_bar_init)
00716     );
00717   // Set initial condition for x_bar_dot = [ x_dot; s_bar_dot ]
00718   state_and_sens_ic.set_x_dot(
00719     stateAndSensModel->create_x_bar_vec(state_ic.get_x_dot(), s_bar_dot_init)
00720     );
00721 
00722   return state_and_sens_ic;
00723 
00724 }
00725 
00726 
00732 template<class Scalar>
00733 inline
00734 Thyra::ModelEvaluatorBase::InArgs<Scalar>
00735 extractStateInitialCondition(
00736   const ForwardSensitivityStepper<Scalar> &fwdSensStepper,
00737   const Thyra::ModelEvaluatorBase::InArgs<Scalar> &state_and_sens_ic
00738   )
00739 {
00740 
00741   using Thyra::productVectorBase;
00742   typedef Thyra::ModelEvaluatorBase MEB;
00743 
00744   MEB::InArgs<Scalar>
00745     state_ic = fwdSensStepper.getStateModel()->createInArgs();
00746 
00747   // Copy time, parameters etc.
00748   state_ic.setArgs(state_and_sens_ic);
00749   state_ic.set_x(
00750     productVectorBase(state_and_sens_ic.get_x())->getVectorBlock(0));
00751   state_ic.set_x_dot(
00752     productVectorBase(state_and_sens_ic.get_x_dot())->getVectorBlock(0));
00753 
00754   return state_ic;
00755 
00756 }
00757 
00758 
00759 //
00760 // Implementation
00761 //
00762 
00763 
00764 // Static members
00765 
00766 
00767 template<class Scalar>
00768 const std::string ForwardSensitivityStepper<Scalar>::forceUpToDateW_name_
00769 = "Force Up-To-Date Jacobian";
00770 
00771 template<class Scalar>
00772 const bool ForwardSensitivityStepper<Scalar>::forceUpToDateW_default_
00773 = true;
00774 
00775 
00776 // Constructors, Intializers, Misc.
00777 
00778 
00779 template<class Scalar>
00780 ForwardSensitivityStepper<Scalar>::ForwardSensitivityStepper()
00781   :forceUpToDateW_(false),
00782    isSingleResidualStepper_(false)
00783 {}
00784 
00785 
00786 template<class Scalar>
00787 void ForwardSensitivityStepper<Scalar>::initializeSyncedSteppers(
00788   const RCP<const Thyra::ModelEvaluator<Scalar> > &stateModel,
00789   const int p_index,
00790   const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateBasePoint,
00791   const RCP<StepperBase<Scalar> > &stateStepper,
00792   const RCP<Thyra::NonlinearSolverBase<Scalar> > &stateTimeStepSolver,
00793   const RCP<StepperBase<Scalar> > &sensStepper,
00794   const RCP<Thyra::NonlinearSolverBase<Scalar> > &sensTimeStepSolver
00795   )
00796 
00797 {
00798   initializeCommon( stateModel, p_index, Teuchos::null, stateBasePoint, stateStepper,
00799     stateTimeStepSolver, sensStepper, sensTimeStepSolver );
00800 }
00801 
00802 
00803 template<class Scalar>
00804 void ForwardSensitivityStepper<Scalar>::initializeSyncedSteppersInitCondOnly(
00805   const RCP<const Thyra::ModelEvaluator<Scalar> >& stateModel,
00806   const RCP<const Thyra::VectorSpaceBase<Scalar> >& p_space,
00807   const Thyra::ModelEvaluatorBase::InArgs<Scalar>& stateBasePoint,
00808   const RCP<StepperBase<Scalar> >& stateStepper,
00809   const RCP<Thyra::NonlinearSolverBase<Scalar> >& stateTimeStepSolver,
00810   const RCP<StepperBase<Scalar> >& sensStepper,
00811   const RCP<Thyra::NonlinearSolverBase<Scalar> >& sensTimeStepSolver
00812   )
00813 {
00814   initializeCommon(stateModel, -1, p_space, stateBasePoint, stateStepper,
00815     stateTimeStepSolver, sensStepper, sensTimeStepSolver );
00816 }
00817 
00818 
00819 template<class Scalar>
00820 void ForwardSensitivityStepper<Scalar>::initializeDecoupledSteppers(
00821   const RCP<const Thyra::ModelEvaluator<Scalar> > &stateModel,
00822   const int p_index,
00823   const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateBasePoint,
00824   const RCP<StepperBase<Scalar> > &stateStepper,
00825   const RCP<Thyra::NonlinearSolverBase<Scalar> > &stateTimeStepSolver,
00826   const RCP<IntegratorBase<Scalar> > &stateIntegrator,
00827   const Scalar &finalTime,
00828   const RCP<StepperBase<Scalar> > &sensStepper,
00829   const RCP<Thyra::NonlinearSolverBase<Scalar> > &sensTimeStepSolver
00830   )
00831 {
00832   TEUCHOS_ASSERT(nonnull(stateIntegrator));
00833   initializeCommon( stateModel, p_index, Teuchos::null, stateBasePoint, stateStepper,
00834     stateTimeStepSolver, sensStepper, sensTimeStepSolver );
00835   stateIntegrator_ = stateIntegrator;
00836   finalTime_ = finalTime;
00837 }
00838 
00839 
00840 template<class Scalar>
00841 bool ForwardSensitivityStepper<Scalar>::stateModelIsConst() const
00842 {
00843   return stateModel_.isConst();
00844 }
00845 
00846 
00847 template<class Scalar>
00848 RCP<const Thyra::ModelEvaluator<Scalar> >
00849 ForwardSensitivityStepper<Scalar>::getStateModel() const
00850 {
00851   return stateModel_.getConstObj();
00852 }
00853 
00854 
00855 template<class Scalar>
00856 RCP<StepperBase<Scalar> >
00857 ForwardSensitivityStepper<Scalar>::getNonconstStateStepper()
00858 {
00859   return stateStepper_;
00860 }
00861 
00862 
00863 template<class Scalar>
00864 RCP<const ForwardSensitivityModelEvaluatorBase<Scalar> >
00865 ForwardSensitivityStepper<Scalar>::getFwdSensModel() const
00866 {
00867   return sensModel_;
00868 }
00869 
00870 
00871 template<class Scalar>
00872 RCP<const StateAndForwardSensitivityModelEvaluator<Scalar> >
00873 ForwardSensitivityStepper<Scalar>::getStateAndFwdSensModel() const
00874 {
00875   return stateAndSensModel_;
00876 }
00877 
00878 
00879 // Overridden from Teuchos::ParameterListAcceptor
00880 
00881 
00882 template<class Scalar>
00883 void ForwardSensitivityStepper<Scalar>::setParameterList(
00884   RCP<Teuchos::ParameterList> const& paramList
00885   )
00886 {
00887   TEUCHOS_TEST_FOR_EXCEPT(is_null(paramList));
00888   paramList->validateParameters(*getValidParameters());
00889   this->setMyParamList(paramList);
00890   forceUpToDateW_ = paramList->get(forceUpToDateW_name_,forceUpToDateW_default_);
00891   Teuchos::readVerboseObjectSublist(&*paramList,this);
00892 }
00893 
00894 
00895 template<class Scalar>
00896 RCP<const Teuchos::ParameterList>
00897 ForwardSensitivityStepper<Scalar>::getValidParameters() const
00898 {
00899   static RCP<const ParameterList> validPL;
00900   if (is_null(validPL) ) {
00901     RCP<ParameterList> pl = Teuchos::parameterList();
00902     pl->set( forceUpToDateW_name_, forceUpToDateW_default_,
00903       "If set to true, then the Jacobian matrix W used in the\n"
00904       "state timestep equation will be forced to be up to date\n"
00905       "with the final value for x for the nonlinear solve.  If\n"
00906       "you are willing to live with slightly less accurate sensitivities\n"
00907       "then set this to false."
00908       );
00909     Teuchos::setupVerboseObjectSublist(&*pl);
00910     validPL = pl;
00911   }
00912   return validPL;
00913 }
00914 
00915 
00916 // Overridden from StepperBase
00917 
00918 template<class Scalar>
00919 bool ForwardSensitivityStepper<Scalar>::acceptsModel() const
00920 {
00921   return false;
00922 }
00923 
00924 template<class Scalar>
00925 void ForwardSensitivityStepper<Scalar>::setModel(
00926   const RCP<const Thyra::ModelEvaluator<Scalar> >& model
00927   )
00928 {
00929   TEUCHOS_TEST_FOR_EXCEPT_MSG( true,
00930     "Error, this stepper subclass does not accept a model"
00931     " as defined by the StepperBase interface!");
00932 }
00933 
00934 
00935 template<class Scalar>
00936 void ForwardSensitivityStepper<Scalar>::setNonconstModel(
00937   const RCP<Thyra::ModelEvaluator<Scalar> >& model
00938   )
00939 {
00940   TEUCHOS_TEST_FOR_EXCEPT_MSG( true,
00941     "Error, this stepper subclass does not accept a model"
00942     " as defined by the StepperBase interface!");
00943 }
00944 
00945 
00946 template<class Scalar>
00947 RCP<const Thyra::ModelEvaluator<Scalar> >
00948 ForwardSensitivityStepper<Scalar>::getModel() const
00949 {
00950   return stateAndSensModel_;
00951 }
00952 
00953 
00954 template<class Scalar>
00955 RCP<Thyra::ModelEvaluator<Scalar> >
00956 ForwardSensitivityStepper<Scalar>::getNonconstModel()
00957 {
00958   return stateAndSensModel_;
00959 }
00960 
00961 
00962 template<class Scalar>
00963 void ForwardSensitivityStepper<Scalar>::setInitialCondition(
00964   const Thyra::ModelEvaluatorBase::InArgs<Scalar> &state_and_sens_ic
00965   )
00966 {
00967 
00968   typedef Thyra::ModelEvaluatorBase MEB;
00969 
00970   stateAndSensBasePoint_ = state_and_sens_ic;
00971 
00972   // Get the product vectors for x_bar = [ x; s_bar ] and x_bar_dot
00973 
00974   TEUCHOS_TEST_FOR_EXCEPTION(
00975     is_null(state_and_sens_ic.get_x()), std::logic_error,
00976     "Error, the initial condition for x_bar = [ x; s_bar ] can not be null!" );
00977 
00978   const RCP<const Thyra::ProductVectorBase<Scalar> >
00979     x_bar_init = Thyra::productVectorBase<Scalar>(
00980       state_and_sens_ic.get_x()
00981       );
00982 
00983   RCP<const Thyra::ProductVectorBase<Scalar> > x_bar_dot_init;
00984   if (state_and_sens_ic.supports(MEB::IN_ARG_x_dot)) {
00985       x_bar_dot_init = Thyra::productVectorBase<Scalar>(
00986         state_and_sens_ic.get_x_dot()
00987         );
00988   }
00989 
00990   // Remove x and x_dot from state_and_sens_ic_in to avoid cloning x and x dot!
00991 
00992   Thyra::ModelEvaluatorBase::InArgs<Scalar>
00993     state_and_sens_ic_no_x = state_and_sens_ic;
00994   state_and_sens_ic_no_x.set_x(Teuchos::null);
00995   if (state_and_sens_ic_no_x.supports(MEB::IN_ARG_x_dot)) {
00996     state_and_sens_ic_no_x.set_x_dot(Teuchos::null);
00997   }
00998 
00999   // Set initial condition for the state
01000 
01001   MEB::InArgs<Scalar> state_ic = stateModel_->createInArgs();
01002   state_ic.setArgs(state_and_sens_ic_no_x,true,true); // Set time, parameters etc.
01003   state_ic.set_x(x_bar_init->getVectorBlock(0)->clone_v());
01004   if (state_ic.supports(MEB::IN_ARG_x_dot)) {
01005     state_ic.set_x_dot(
01006         !is_null(x_bar_dot_init)
01007         ? x_bar_dot_init->getVectorBlock(0)->clone_v()
01008         : Teuchos::null
01009         );
01010   }
01011   stateStepper_->setInitialCondition(state_ic);
01012 
01013   // Set up the integrator if needed
01014   //if (!is_null(stateIntegrator_)) {
01015   //  stateIntegrator_->setStepper( stateStepper_, finalTime_ );
01016   //  sensModel_->setStateIntegrator( stateIntegrator_, state_ic );
01017   //}
01018 
01019   // Set initial condition for the sensitivities
01020 
01021   MEB::InArgs<Scalar> sens_ic = sensModel_->createInArgs();
01022   sens_ic.setArgs(state_and_sens_ic_no_x,true,true); // Set time etc.
01023   sens_ic.set_x(x_bar_init->getVectorBlock(1)->clone_v());
01024   if (sens_ic.supports(MEB::IN_ARG_x_dot)) {
01025     sens_ic.set_x_dot(
01026         !is_null(x_bar_dot_init)
01027         ? x_bar_dot_init->getVectorBlock(1)->clone_v()
01028         : Teuchos::null
01029         );
01030   }
01031   sensStepper_->setInitialCondition(sens_ic);
01032 
01033 }
01034 
01035 
01036 template<class Scalar>
01037 Thyra::ModelEvaluatorBase::InArgs<Scalar>
01038 ForwardSensitivityStepper<Scalar>::getInitialCondition() const
01039 {
01040   return stateAndSensBasePoint_;
01041 }
01042 
01043 
01044 template<class Scalar>
01045 Scalar
01046 ForwardSensitivityStepper<Scalar>::takeStep(
01047   Scalar dt, StepSizeType stepType
01048   )
01049 {
01050 
01051   RYTHMOS_FUNC_TIME_MONITOR("Rythmos:ForwardSensitivityStepper::takeStep");
01052 
01053   if (!is_null(stateIntegrator_)) {
01054     return takeDecoupledStep(dt,stepType);
01055   }
01056 
01057   return takeSyncedStep(dt,stepType);
01058 
01059 }
01060 
01061 
01062 template<class Scalar>
01063 const StepStatus<Scalar>
01064 ForwardSensitivityStepper<Scalar>::getStepStatus() const
01065 {
01066 
01067   const StepStatus<Scalar> sensStepStatus = sensStepper_->getStepStatus();
01068   StepStatus<Scalar> stepStatus;
01069 
01070   stepStatus.message = sensStepStatus.message;
01071   stepStatus.stepStatus = sensStepStatus.stepStatus;
01072   stepStatus.stepLETStatus = sensStepStatus.stepLETStatus;
01073   stepStatus.stepSize = sensStepStatus.stepSize;
01074   stepStatus.order = sensStepStatus.order;
01075   stepStatus.time = sensStepStatus.time;
01076   stepStatus.stepLETValue = sensStepStatus.stepLETValue;
01077   stepStatus.extraParameters = sensStepStatus.extraParameters;
01078 
01079   if (is_null(stateIntegrator_)) {
01080     const StepStatus<Scalar>
01081       stateStepStatus = stateStepper_->getStepStatus();
01082     if (!is_null(stateStepStatus.solution) && !is_null(sensStepStatus.solution))
01083       stepStatus.solution = stateAndSensModel_->create_x_bar_vec(
01084         stateStepStatus.solution, sensStepStatus.solution
01085         );
01086     if (!is_null(stateStepStatus.solutionDot) && !is_null(sensStepStatus.solutionDot))
01087       stepStatus.solutionDot = stateAndSensModel_->create_x_bar_vec(
01088         stateStepStatus.solutionDot, sensStepStatus.solutionDot
01089         );
01090   }
01091 
01092   return stepStatus;
01093 
01094 }
01095 
01096 
01097 // Overridden from InterpolationBufferBase
01098 
01099 
01100 template<class Scalar>
01101 RCP<const Thyra::VectorSpaceBase<Scalar> >
01102 ForwardSensitivityStepper<Scalar>::get_x_space() const
01103 {
01104   return stateAndSensModel_->get_x_space();
01105 }
01106 
01107 
01108 template<class Scalar>
01109 void ForwardSensitivityStepper<Scalar>::addPoints(
01110   const Array<Scalar>& time_vec,
01111   const Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >& x_vec,
01112   const Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >& xdot_vec
01113   )
01114 {
01115   TEUCHOS_TEST_FOR_EXCEPT("Not implemented addPoints(...) yet but we could if we wanted!");
01116 }
01117 
01118 
01119 template<class Scalar>
01120 TimeRange<Scalar>
01121 ForwardSensitivityStepper<Scalar>::getTimeRange() const
01122 {
01123   return sensStepper_->getTimeRange();
01124 }
01125 
01126 
01127 template<class Scalar>
01128 void ForwardSensitivityStepper<Scalar>::getPoints(
01129   const Array<Scalar>& time_vec,
01130   Array<RCP<const Thyra::VectorBase<Scalar> > >* x_bar_vec,
01131   Array<RCP<const Thyra::VectorBase<Scalar> > >* x_bar_dot_vec,
01132   Array<ScalarMag>* accuracy_vec
01133   ) const
01134 {
01135 
01136   using Teuchos::as;
01137 
01138 #ifdef HAVE_RYTHMOS_DEBUG
01139   TEUCHOS_TEST_FOR_EXCEPT( as<int>(time_vec.size()) == 0 );
01140 #endif
01141 
01142   const int numTimePoints = time_vec.size();
01143 
01144   if (x_bar_vec)
01145     x_bar_vec->clear();
01146 
01147   if (x_bar_dot_vec)
01148     x_bar_dot_vec->clear();
01149 
01150   Array<RCP<const Thyra::VectorBase<Scalar> > >
01151     x_vec, x_dot_vec;
01152 
01153   if (!is_null(stateIntegrator_)) {
01154     stateIntegrator_->getPoints(
01155       time_vec,
01156       x_bar_vec ? &x_vec: 0,
01157       x_bar_dot_vec ? &x_dot_vec: 0,
01158       0 // Ignoring accuracy from state for now!
01159       );
01160   }
01161   else {
01162     stateStepper_->getPoints(
01163       time_vec,
01164       x_bar_vec ? &x_vec: 0,
01165       x_bar_dot_vec ? &x_dot_vec: 0,
01166       0 // Ignoring accuracy from state for now!
01167       );
01168   }
01169 
01170   Array<RCP<const Thyra::VectorBase<Scalar> > >
01171     s_bar_vec, s_bar_dot_vec;
01172 
01173   sensStepper_->getPoints(
01174     time_vec,
01175     x_bar_vec ? &s_bar_vec: 0,
01176     x_bar_dot_vec ? &s_bar_dot_vec: 0,
01177     accuracy_vec
01178     );
01179 
01180   if ( x_bar_vec ) {
01181     for ( int i = 0; i < numTimePoints; ++i ) {
01182       x_bar_vec->push_back(
01183         stateAndSensModel_->create_x_bar_vec(x_vec[i],s_bar_vec[i])
01184         );
01185     }
01186   }
01187 
01188   if ( x_bar_dot_vec ) {
01189     for ( int i = 0; i < numTimePoints; ++i ) {
01190       x_bar_dot_vec->push_back(
01191         stateAndSensModel_->create_x_bar_vec(x_dot_vec[i],s_bar_dot_vec[i])
01192         );
01193     }
01194   }
01195 
01196 }
01197 
01198 
01199 template<class Scalar>
01200 void ForwardSensitivityStepper<Scalar>::getNodes(
01201   Array<Scalar>* time_vec
01202   ) const
01203 {
01204   TEUCHOS_ASSERT( time_vec != NULL );
01205   time_vec->clear();
01206   if (is_null(stateIntegrator_) && is_null(stateStepper_)) {
01207     return;
01208   }
01209   if (!is_null(stateIntegrator_)) {
01210     stateIntegrator_->getNodes(time_vec);
01211   }
01212   else {
01213     stateStepper_->getNodes(time_vec);
01214   }
01215 }
01216 
01217 
01218 template<class Scalar>
01219 void ForwardSensitivityStepper<Scalar>::removeNodes(
01220   Array<Scalar>& time_vec
01221   )
01222 {
01223   TEUCHOS_TEST_FOR_EXCEPT("Not implemented yet but we can!");
01224 }
01225 
01226 
01227 template<class Scalar>
01228 int ForwardSensitivityStepper<Scalar>::getOrder() const
01229 {
01230   return sensStepper_->getOrder();
01231   // Note: This assumes that stateStepper will have the same order!
01232 }
01233 
01234 
01235 // private
01236 
01237 
01238 template<class Scalar>
01239 void ForwardSensitivityStepper<Scalar>::initializeCommon(
01240   const RCP<const Thyra::ModelEvaluator<Scalar> >& stateModel,
01241   const int p_index,
01242   const RCP<const Thyra::VectorSpaceBase<Scalar> > &p_space,
01243   const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateBasePoint,
01244   const RCP<StepperBase<Scalar> > &stateStepper,
01245   const RCP<Thyra::NonlinearSolverBase<Scalar> > &stateTimeStepSolver,
01246   const RCP<StepperBase<Scalar> > &sensStepper,
01247   const RCP<Thyra::NonlinearSolverBase<Scalar> > &sensTimeStepSolver
01248   )
01249 {
01250 
01251   using Teuchos::rcp_implicit_cast;
01252   using Teuchos::rcp_dynamic_cast;
01253 
01254   typedef Thyra::ModelEvaluatorBase MEB;
01255 
01256   //
01257   // Validate input
01258   //
01259 
01260   TEUCHOS_ASSERT( p_index >= 0 || nonnull(p_space) );
01261   if (nonnull(p_space)) {
01262     TEUCHOS_ASSERT_EQUALITY(p_index, -1);
01263   }
01264   if (p_index >= 0) {
01265     TEUCHOS_ASSERT(is_null(p_space));
01266   }
01267   TEUCHOS_TEST_FOR_EXCEPT( is_null(stateModel) );
01268   TEUCHOS_TEST_FOR_EXCEPT( is_null(stateStepper) );
01269   if (stateStepper->isImplicit()) {
01270     TEUCHOS_TEST_FOR_EXCEPT( is_null(stateTimeStepSolver) ); // allow to be null for explicit methods
01271   }
01272 
01273   //
01274   // Create the sensModel which will do some more validation
01275   //
01276 
01277   RCP<ForwardSensitivityModelEvaluatorBase<Scalar> > sensModel;
01278   MEB::InArgs<Scalar> stateModelInArgs = stateModel->createInArgs();
01279   if (stateModelInArgs.supports(MEB::IN_ARG_x_dot)) {
01280     // Implicit DE formulation
01281     sensModel = Teuchos::rcp(new ForwardSensitivityImplicitModelEvaluator<Scalar>);
01282   }
01283   else {
01284     // Explicit DE formulation
01285     sensModel = Teuchos::rcp(new ForwardSensitivityExplicitModelEvaluator<Scalar>);
01286   }
01287 
01288   if (p_index >= 0) {
01289     sensModel->initializeStructure(stateModel, p_index);
01290   }
01291   else {
01292     sensModel->initializeStructureInitCondOnly(stateModel, p_space);
01293   }
01294 
01295   //
01296   // Get the input objects
01297   //
01298 
01299   stateModel_.initialize(stateModel);
01300 
01301   stateBasePoint_ = stateBasePoint;
01302 
01303   stateStepper_ = stateStepper;
01304 
01305   stateTimeStepSolver_ = stateTimeStepSolver;
01306 
01307   sensModel_ = sensModel;
01308 
01309   stateAndSensModel_ = Teuchos::rcp(new StateAndForwardSensitivityModelEvaluator<Scalar>);
01310   stateAndSensModel_->initializeStructure(sensModel_);
01311 
01312   if (!is_null(sensStepper)) {
01313     sensStepper_ = sensStepper;
01314   }
01315   else {
01316     sensStepper_ = stateStepper_->cloneStepperAlgorithm();
01317     TEUCHOS_TEST_FOR_EXCEPTION(
01318       is_null(sensStepper_), std::logic_error,
01319       "Error, if the client does not pass in a stepper for the senitivity\n"
01320       "equations then the stateStepper object must support cloning to create\n"
01321       "the sensitivity stepper!"
01322       );
01323   }
01324 
01325   if (!is_null(sensTimeStepSolver)) {
01326     sensTimeStepSolver_ = sensTimeStepSolver;
01327   }
01328   else {
01329     RCP<Thyra::LinearNonlinearSolver<Scalar> >
01330       linearNonlinearSolver(new Thyra::LinearNonlinearSolver<Scalar>);
01331     // ToDo: Set tolerance on the nonlinear solver???
01332     sensTimeStepSolver_ = linearNonlinearSolver;
01333   }
01334 
01335   //
01336   // Setup the steppers
01337   //
01338 
01339   isSingleResidualStepper_ = true; // ToDo: Add dynamic cast on
01340                                    // stateTimeStepSolver to check this!
01341 
01342   setStepperModel(Teuchos::inOutArg(*stateStepper_),stateModel_);
01343   if (stateStepper_->isImplicit()) {
01344     rcp_dynamic_cast<SolverAcceptingStepperBase<Scalar> >(
01345         stateStepper_,true)->setSolver(stateTimeStepSolver_);
01346   }
01347   sensStepper_->setModel(sensModel_);
01348   if (sensStepper_->isImplicit()) {
01349     rcp_dynamic_cast<SolverAcceptingStepperBase<Scalar> >(
01350         sensStepper_,true)->setSolver(sensTimeStepSolver_);
01351   }
01352 
01353   stateBasePoint_t_ = stateModel_->createInArgs();
01354 
01355   // 2007/05/18: rabartl: ToDo: Move the above initialization code to give
01356   // setInitializeCondition(...) a chance to set the initial condition.
01357 
01358 }
01359 
01360 
01361 template<class Scalar>
01362 Scalar ForwardSensitivityStepper<Scalar>::takeSyncedStep(
01363   Scalar dt, StepSizeType stepType
01364   )
01365 {
01366 
01367   RYTHMOS_FUNC_TIME_MONITOR_DIFF("Rythmos:ForwardSensitivityStepper::takeStep: synced",
01368     TopLevel);
01369 
01370   using Teuchos::as;
01371   typedef Teuchos::ScalarTraits<Scalar> ST;
01372   typedef Teuchos::VerboseObjectTempState<InterpolationBufferBase<Scalar> > VOTSIBB;
01373   typedef Thyra::ModelEvaluatorBase MEB;
01374 
01375   RCP<Teuchos::FancyOStream> out = this->getOStream();
01376   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
01377   const bool lowTrace =
01378     ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) );
01379   const bool mediumTrace =
01380     ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_MEDIUM) );
01381   Teuchos::OSTab tab(out);
01382 
01383   if (lowTrace) {
01384     *out
01385       << "\nEntering " << TypeNameTraits<ForwardSensitivityStepper<Scalar> >::name()
01386       << "::takeSyncedStep("<<dt<<","<<toString(stepType)<<") ...\n";
01387   }
01388 
01389   //
01390   // A) Compute the state timestep
01391   //
01392 
01393   if (lowTrace) {
01394     *out
01395       << "\nTaking state step using stepper : "
01396       << stateStepper_->description() << "\n";
01397   }
01398 
01399   Scalar state_dt = -1.0;
01400   {
01401     RYTHMOS_FUNC_TIME_MONITOR("Rythmos:ForwardSensitivityStepper::takeStep: stateStep");
01402     VOTSIBB stateStepper_outputTempState(stateStepper_,out,verbLevel);
01403     state_dt = stateStepper_->takeStep(dt,stepType);
01404   }
01405 
01406   if (state_dt < Scalar(-ST::one())) {
01407     if (lowTrace)
01408       *out << "\nThe state stepper has failed so return a failed timestep!\n";
01409     return state_dt;
01410   }
01411 
01412   {
01413     RYTHMOS_FUNC_TIME_MONITOR("Rythmos:ForwardSensitivityStepper::takeStep: updateSensModel");
01414     // Set up the sensitivity model for this timestep
01415     sensModel_->initializePointState(Teuchos::inOutArg(*stateStepper_),forceUpToDateW_);
01416   }
01417 
01418   //
01419   // C) Compute the sensitivity timestep for the exact same timestep as was
01420   // used for the state solve.
01421   //
01422 
01423   if (lowTrace) {
01424     *out
01425       << "\nTaking sensitivity step using stepper : "
01426       << sensStepper_->description() << "\n";
01427   }
01428 
01429   Scalar sens_dt = -1.0;
01430   {
01431     RYTHMOS_FUNC_TIME_MONITOR("Rythmos:ForwardSensitivityStepper::takeStep: sensStep");
01432     // Copy the step control data to make sure that the sensStepper takes the
01433     // same type of step that the statStepper took.  This is needed to ensure
01434     // that the W matrix is the same for one.
01435     sensStepper_->setStepControlData(*stateStepper_);
01436     VOTSIBB sensStepper_outputTempState(sensStepper_,out,verbLevel);
01437     sens_dt = sensStepper_->takeStep(state_dt,STEP_TYPE_FIXED);
01438   }
01439 
01440   if (mediumTrace) {
01441     const StepStatus<Scalar> sensStepStatus = sensStepper_->getStepStatus();
01442     *out << "\nSensitivity step status:\n" << sensStepStatus;
01443   }
01444 
01445   TEUCHOS_TEST_FOR_EXCEPTION(
01446     sens_dt != state_dt, std::logic_error,
01447     "Error, the sensitivity step failed for some reason.  We should\n"
01448     "just return a negative step size and reject the step but currently\n"
01449     "there is no way to roll back the state timestep it for back to\n"
01450     "the status before this function was called!"
01451     );
01452 
01453   // 2007/05/18: rabartl: ToDo: If stepType == STEP_TYPE_VARIABLE and the state
01454   // timestep sucessed but the sensitivity timestep failed, then we need to
01455   // really throw an excpetion because there is nothing that we can really do
01456   // here!
01457 
01458   // 2007/05/18: rabartl: ToDo: Replace the above std::logic_error type with
01459   // a Rythmos::CatastrophicFailure or just use Thyra::CatastrophicFailure!
01460 
01461   if (lowTrace) {
01462     *out
01463       << "\nLeaving " << TypeNameTraits<ForwardSensitivityStepper<Scalar> >::name()
01464       << "::takeSyncedStep("<<dt<<","<<toString(stepType)<<") ...\n";
01465   }
01466 
01467   return state_dt;
01468 
01469 }
01470 
01471 
01472 template<class Scalar>
01473 Scalar ForwardSensitivityStepper<Scalar>::takeDecoupledStep(
01474   Scalar dt, StepSizeType stepType
01475   )
01476 {
01477 
01478   RYTHMOS_FUNC_TIME_MONITOR("Rythmos:ForwardSensitivityStepper::takeStep: decoupled");
01479 
01480   using Teuchos::as;
01481   typedef Teuchos::ScalarTraits<Scalar> ST;
01482   typedef Teuchos::VerboseObjectTempState<InterpolationBufferBase<Scalar> > VOTSIBB;
01483   typedef Thyra::ModelEvaluatorBase MEB;
01484 
01485   RCP<Teuchos::FancyOStream> out = this->getOStream();
01486   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
01487   const bool lowTrace =
01488     ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) );
01489   const bool mediumTrace =
01490     ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_MEDIUM) );
01491   Teuchos::OSTab tab(out);
01492 
01493   if (lowTrace) {
01494     *out
01495       << "\nEntering " << TypeNameTraits<ForwardSensitivityStepper<Scalar> >::name()
01496       << "::takeDecoupledStep("<<dt<<","<<toString(stepType)<<") ...\n";
01497   }
01498 
01499   //
01500   // A) Take the sens timestep
01501   //
01502 
01503   if (lowTrace) {
01504     *out
01505       << "\nTaking sensitivity step using stepper : "
01506       << sensStepper_->description() << "\n";
01507   }
01508 
01509   Scalar sens_dt = -1.0;
01510   VOTSIBB sensStepper_outputTempState(sensStepper_,out,verbLevel);
01511   sens_dt = sensStepper_->takeStep(dt,stepType);
01512 
01513   if (mediumTrace) {
01514     const StepStatus<Scalar> sensStepStatus = sensStepper_->getStepStatus();
01515     *out << "\nSensitivity step status:\n" << sensStepStatus;
01516   }
01517 
01518   //
01519   // B) Wipe out all state interp buffer info before this sens timestep
01520   //
01521 
01522   //TEUCHOS_TEST_FOR_EXCEPT(true);
01523 
01524   if (lowTrace) {
01525     *out
01526       << "\nLeaving " << TypeNameTraits<ForwardSensitivityStepper<Scalar> >::name()
01527       << "::takeDecoupledStep("<<dt<<","<<toString(stepType)<<") ...\n";
01528   }
01529 
01530   return sens_dt;
01531 
01532 }
01533 
01534 
01535 } // namespace Rythmos
01536 
01537 
01538 #endif //RYTHMOS_FORWARD_SENSITIVITY_STEPPER_HPP
 All Classes Functions Variables Typedefs Friends