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   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   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   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   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 #ifdef ENABLE_RYTHMOS_TIMERS
01052   TEUCHOS_FUNC_TIME_MONITOR("Rythmos:ForwardSensitivityStepper::takeStep");
01053 #endif
01054 
01055   if (!is_null(stateIntegrator_)) {
01056     return takeDecoupledStep(dt,stepType);
01057   }
01058 
01059   return takeSyncedStep(dt,stepType);
01060 
01061 }
01062 
01063 
01064 template<class Scalar> 
01065 const StepStatus<Scalar>
01066 ForwardSensitivityStepper<Scalar>::getStepStatus() const
01067 {
01068 
01069   const StepStatus<Scalar> sensStepStatus = sensStepper_->getStepStatus();
01070   StepStatus<Scalar> stepStatus;
01071   
01072   stepStatus.message = sensStepStatus.message;
01073   stepStatus.stepStatus = sensStepStatus.stepStatus;
01074   stepStatus.stepLETStatus = sensStepStatus.stepLETStatus;
01075   stepStatus.stepSize = sensStepStatus.stepSize;
01076   stepStatus.order = sensStepStatus.order;
01077   stepStatus.time = sensStepStatus.time;
01078   stepStatus.stepLETValue = sensStepStatus.stepLETValue;
01079   stepStatus.extraParameters = sensStepStatus.extraParameters;
01080 
01081   if (is_null(stateIntegrator_)) {
01082     const StepStatus<Scalar> 
01083       stateStepStatus = stateStepper_->getStepStatus();
01084     if (!is_null(stateStepStatus.solution) && !is_null(sensStepStatus.solution))
01085       stepStatus.solution = stateAndSensModel_->create_x_bar_vec(
01086         stateStepStatus.solution, sensStepStatus.solution
01087         );
01088     if (!is_null(stateStepStatus.solutionDot) && !is_null(sensStepStatus.solutionDot))
01089       stepStatus.solutionDot = stateAndSensModel_->create_x_bar_vec(
01090         stateStepStatus.solutionDot, sensStepStatus.solutionDot
01091         );
01092   }
01093 
01094   return stepStatus;
01095 
01096 }
01097 
01098 
01099 // Overridden from InterpolationBufferBase
01100 
01101 
01102 template<class Scalar> 
01103 RCP<const Thyra::VectorSpaceBase<Scalar> >
01104 ForwardSensitivityStepper<Scalar>::get_x_space() const
01105 {
01106   return stateAndSensModel_->get_x_space();
01107 }
01108 
01109 
01110 template<class Scalar> 
01111 void ForwardSensitivityStepper<Scalar>::addPoints(
01112   const Array<Scalar>& time_vec,
01113   const Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >& x_vec,
01114   const Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >& xdot_vec
01115   )
01116 {
01117   TEST_FOR_EXCEPT("Not implemented addPoints(...) yet but we could if we wanted!");
01118 }
01119 
01120 
01121 template<class Scalar> 
01122 TimeRange<Scalar>
01123 ForwardSensitivityStepper<Scalar>::getTimeRange() const
01124 {
01125   return sensStepper_->getTimeRange();
01126 }
01127 
01128 
01129 template<class Scalar> 
01130 void ForwardSensitivityStepper<Scalar>::getPoints(
01131   const Array<Scalar>& time_vec,
01132   Array<RCP<const Thyra::VectorBase<Scalar> > >* x_bar_vec,
01133   Array<RCP<const Thyra::VectorBase<Scalar> > >* x_bar_dot_vec,
01134   Array<ScalarMag>* accuracy_vec
01135   ) const
01136 {
01137 
01138   using Teuchos::as;
01139 
01140 #ifdef RYTHMOS_DEBUG
01141   TEST_FOR_EXCEPT( as<int>(time_vec.size()) == 0 );
01142 #endif
01143 
01144   const int numTimePoints = time_vec.size();
01145 
01146   if (x_bar_vec)
01147     x_bar_vec->clear();
01148 
01149   if (x_bar_dot_vec)
01150     x_bar_dot_vec->clear();
01151   
01152   Array<RCP<const Thyra::VectorBase<Scalar> > >
01153     x_vec, x_dot_vec;
01154 
01155   if (!is_null(stateIntegrator_)) {
01156     stateIntegrator_->getPoints(
01157       time_vec,
01158       x_bar_vec ? &x_vec: 0,
01159       x_bar_dot_vec ? &x_dot_vec: 0,
01160       0 // Ignoring accuracy from state for now!
01161       );
01162   }
01163   else {
01164     stateStepper_->getPoints(
01165       time_vec,
01166       x_bar_vec ? &x_vec: 0,
01167       x_bar_dot_vec ? &x_dot_vec: 0,
01168       0 // Ignoring accuracy from state for now!
01169       );
01170   }
01171   
01172   Array<RCP<const Thyra::VectorBase<Scalar> > >
01173     s_bar_vec, s_bar_dot_vec;
01174 
01175   sensStepper_->getPoints(
01176     time_vec,
01177     x_bar_vec ? &s_bar_vec: 0,
01178     x_bar_dot_vec ? &s_bar_dot_vec: 0,
01179     accuracy_vec
01180     );
01181   
01182   if ( x_bar_vec ) {
01183     for ( int i = 0; i < numTimePoints; ++i ) {
01184       x_bar_vec->push_back(
01185         stateAndSensModel_->create_x_bar_vec(x_vec[i],s_bar_vec[i])
01186         );
01187     }
01188   }
01189   
01190   if ( x_bar_dot_vec ) {
01191     for ( int i = 0; i < numTimePoints; ++i ) {
01192       x_bar_dot_vec->push_back(
01193         stateAndSensModel_->create_x_bar_vec(x_dot_vec[i],s_bar_dot_vec[i])
01194         );
01195     }
01196   }
01197 
01198 }
01199 
01200 
01201 template<class Scalar>
01202 void ForwardSensitivityStepper<Scalar>::getNodes(
01203   Array<Scalar>* time_vec
01204   ) const
01205 {
01206   TEUCHOS_ASSERT( time_vec != NULL );
01207   time_vec->clear();
01208   if (is_null(stateIntegrator_) && is_null(stateStepper_)) {
01209     return;
01210   }
01211   if (!is_null(stateIntegrator_)) {
01212     stateIntegrator_->getNodes(time_vec);
01213   }
01214   else {
01215     stateStepper_->getNodes(time_vec);
01216   }
01217 }
01218 
01219 
01220 template<class Scalar> 
01221 void ForwardSensitivityStepper<Scalar>::removeNodes(
01222   Array<Scalar>& time_vec
01223   )
01224 {
01225   TEST_FOR_EXCEPT("Not implemented yet but we can!");
01226 }
01227 
01228 
01229 template<class Scalar> 
01230 int ForwardSensitivityStepper<Scalar>::getOrder() const
01231 {
01232   return sensStepper_->getOrder();
01233   // Note: This assumes that stateStepper will have the same order!
01234 }
01235 
01236 
01237 // private
01238 
01239 
01240 template<class Scalar> 
01241 void ForwardSensitivityStepper<Scalar>::initializeCommon(
01242   const RCP<const Thyra::ModelEvaluator<Scalar> >& stateModel,
01243   const int p_index,
01244   const RCP<const Thyra::VectorSpaceBase<Scalar> > &p_space,
01245   const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateBasePoint,
01246   const RCP<StepperBase<Scalar> > &stateStepper,
01247   const RCP<Thyra::NonlinearSolverBase<Scalar> > &stateTimeStepSolver,
01248   const RCP<StepperBase<Scalar> > &sensStepper,
01249   const RCP<Thyra::NonlinearSolverBase<Scalar> > &sensTimeStepSolver
01250   )
01251 {
01252 
01253   using Teuchos::rcp_implicit_cast;
01254   using Teuchos::rcp_dynamic_cast;
01255 
01256   typedef Thyra::ModelEvaluatorBase MEB;
01257 
01258   //
01259   // Validate input
01260   //
01261 
01262   TEUCHOS_ASSERT( p_index >= 0 || nonnull(p_space) );
01263   if (nonnull(p_space)) {
01264     TEUCHOS_ASSERT_EQUALITY(p_index, -1);
01265   }
01266   if (p_index >= 0) {
01267     TEUCHOS_ASSERT(is_null(p_space));
01268   }
01269   TEST_FOR_EXCEPT( is_null(stateModel) );
01270   TEST_FOR_EXCEPT( is_null(stateStepper) );
01271   if (stateStepper->isImplicit()) {
01272     TEST_FOR_EXCEPT( is_null(stateTimeStepSolver) ); // allow to be null for explicit methods
01273   }
01274 
01275   //
01276   // Create the sensModel which will do some more validation
01277   //
01278   
01279   RCP<ForwardSensitivityModelEvaluatorBase<Scalar> > sensModel;
01280   MEB::InArgs<Scalar> stateModelInArgs = stateModel->createInArgs();
01281   if (stateModelInArgs.supports(MEB::IN_ARG_x_dot)) {
01282     // Implicit DE formulation
01283     sensModel = Teuchos::rcp(new ForwardSensitivityImplicitModelEvaluator<Scalar>);
01284   }
01285   else {
01286     // Explicit DE formulation
01287     sensModel = Teuchos::rcp(new ForwardSensitivityExplicitModelEvaluator<Scalar>);
01288   }
01289 
01290   if (p_index >= 0) {
01291     sensModel->initializeStructure(stateModel, p_index);
01292   }
01293   else {
01294     sensModel->initializeStructureInitCondOnly(stateModel, p_space);
01295   }
01296   
01297   //
01298   // Get the input objects
01299   //
01300 
01301   stateModel_.initialize(stateModel);
01302 
01303   stateBasePoint_ = stateBasePoint;
01304 
01305   stateStepper_ = stateStepper;
01306   
01307   stateTimeStepSolver_ = stateTimeStepSolver;
01308 
01309   sensModel_ = sensModel;
01310 
01311   stateAndSensModel_ = Teuchos::rcp(new StateAndForwardSensitivityModelEvaluator<Scalar>);
01312   stateAndSensModel_->initializeStructure(sensModel_);
01313 
01314   if (!is_null(sensStepper)) {
01315     sensStepper_ = sensStepper;
01316   }
01317   else {
01318     sensStepper_ = stateStepper_->cloneStepperAlgorithm();
01319     TEST_FOR_EXCEPTION(
01320       is_null(sensStepper_), std::logic_error,
01321       "Error, if the client does not pass in a stepper for the senitivity\n"
01322       "equations then the stateStepper object must support cloning to create\n"
01323       "the sensitivity stepper!"
01324       );
01325   }
01326 
01327   if (!is_null(sensTimeStepSolver)) {
01328     sensTimeStepSolver_ = sensTimeStepSolver;
01329   }
01330   else {
01331     RCP<Thyra::LinearNonlinearSolver<Scalar> >
01332       linearNonlinearSolver(new Thyra::LinearNonlinearSolver<Scalar>);
01333     // ToDo: Set tolerance on the nonlinear solver???
01334     sensTimeStepSolver_ = linearNonlinearSolver;
01335   }
01336 
01337   //
01338   // Setup the steppers
01339   //
01340 
01341   isSingleResidualStepper_ = true; // ToDo: Add dynamic cast on
01342                                    // stateTimeStepSolver to check this!
01343 
01344   setStepperModel(Teuchos::inOutArg(*stateStepper_),stateModel_);
01345   if (stateStepper_->isImplicit()) {
01346     rcp_dynamic_cast<SolverAcceptingStepperBase<Scalar> >(
01347         stateStepper_,true)->setSolver(stateTimeStepSolver_);
01348   }
01349   sensStepper_->setModel(sensModel_); 
01350   if (sensStepper_->isImplicit()) {
01351     rcp_dynamic_cast<SolverAcceptingStepperBase<Scalar> >(
01352         sensStepper_,true)->setSolver(sensTimeStepSolver_);
01353   }
01354 
01355   stateBasePoint_t_ = stateModel_->createInArgs();
01356 
01357   // 2007/05/18: rabartl: ToDo: Move the above initialization code to give
01358   // setInitializeCondition(...) a chance to set the initial condition.
01359 
01360 }
01361 
01362 
01363 template<class Scalar> 
01364 Scalar ForwardSensitivityStepper<Scalar>::takeSyncedStep(
01365   Scalar dt, StepSizeType stepType
01366   )
01367 {
01368 
01369 #ifdef ENABLE_RYTHMOS_TIMERS
01370   TEUCHOS_FUNC_TIME_MONITOR_DIFF("Rythmos:ForwardSensitivityStepper::takeStep: synced",
01371     TopLevel);
01372 #endif
01373 
01374   using Teuchos::as;
01375   typedef Teuchos::ScalarTraits<Scalar> ST;
01376   typedef Teuchos::VerboseObjectTempState<InterpolationBufferBase<Scalar> > VOTSIBB;
01377   typedef Thyra::ModelEvaluatorBase MEB;
01378 
01379   RCP<Teuchos::FancyOStream> out = this->getOStream();
01380   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
01381   const bool lowTrace =
01382     ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) );
01383   const bool mediumTrace =
01384     ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_MEDIUM) );
01385   Teuchos::OSTab tab(out);
01386 
01387   if (lowTrace) {
01388     *out
01389       << "\nEntering " << TypeNameTraits<ForwardSensitivityStepper<Scalar> >::name()
01390       << "::takeSyncedStep("<<dt<<","<<toString(stepType)<<") ...\n"; 
01391   }
01392 
01393   //
01394   // A) Compute the state timestep
01395   //
01396 
01397   if (lowTrace) {
01398     *out
01399       << "\nTaking state step using stepper : "
01400       << stateStepper_->description() << "\n";
01401   }
01402 
01403   Scalar state_dt = -1.0;
01404   {
01405 #ifdef ENABLE_RYTHMOS_TIMERS
01406     TEUCHOS_FUNC_TIME_MONITOR("Rythmos:ForwardSensitivityStepper::takeStep: stateStep");
01407 #endif
01408     VOTSIBB stateStepper_outputTempState(stateStepper_,out,verbLevel);
01409     state_dt = stateStepper_->takeStep(dt,stepType);
01410   }
01411 
01412   if (state_dt < Scalar(-ST::one())) {
01413     if (lowTrace)
01414       *out << "\nThe state stepper has failed so return a failed timestep!\n";
01415     return state_dt;
01416   }
01417 
01418   {
01419 #ifdef ENABLE_RYTHMOS_TIMERS
01420     TEUCHOS_FUNC_TIME_MONITOR("Rythmos:ForwardSensitivityStepper::takeStep: updateSensModel");
01421 #endif
01422     // Set up the sensitivity model for this timestep
01423     sensModel_->initializePointState(Teuchos::inOutArg(*stateStepper_),forceUpToDateW_);
01424   } 
01425 
01426   //
01427   // C) Compute the sensitivity timestep for the exact same timestep as was
01428   // used for the state solve.
01429   //
01430 
01431   if (lowTrace) {
01432     *out
01433       << "\nTaking sensitivity step using stepper : "
01434       << sensStepper_->description() << "\n";
01435   }
01436 
01437   Scalar sens_dt = -1.0;
01438   {
01439 #ifdef ENABLE_RYTHMOS_TIMERS
01440     TEUCHOS_FUNC_TIME_MONITOR("Rythmos:ForwardSensitivityStepper::takeStep: sensStep");
01441 #endif
01442     // Copy the step control data to make sure that the sensStepper takes the
01443     // same type of step that the statStepper took.  This is needed to ensure
01444     // that the W matrix is the same for one.
01445     sensStepper_->setStepControlData(*stateStepper_);
01446     VOTSIBB sensStepper_outputTempState(sensStepper_,out,verbLevel);
01447     sens_dt = sensStepper_->takeStep(state_dt,STEP_TYPE_FIXED);
01448   }
01449 
01450   if (mediumTrace) {
01451     const StepStatus<Scalar> sensStepStatus = sensStepper_->getStepStatus();
01452     *out << "\nSensitivity step status:\n" << sensStepStatus;
01453   }
01454   
01455   TEST_FOR_EXCEPTION(
01456     sens_dt != state_dt, std::logic_error,
01457     "Error, the sensitivity step failed for some reason.  We should\n"
01458     "just return a negative step size and reject the step but currently\n"
01459     "there is no way to roll back the state timestep it for back to\n"
01460     "the status before this function was called!"
01461     );
01462 
01463   // 2007/05/18: rabartl: ToDo: If stepType == STEP_TYPE_VARIABLE and the state
01464   // timestep sucessed but the sensitivity timestep failed, then we need to
01465   // really throw an excpetion because there is nothing that we can really do
01466   // here!
01467 
01468   // 2007/05/18: rabartl: ToDo: Replace the above std::logic_error type with
01469   // a Rythmos::CatastrophicFailure or just use Thyra::CatastrophicFailure!
01470 
01471   if (lowTrace) {
01472     *out
01473       << "\nLeaving " << TypeNameTraits<ForwardSensitivityStepper<Scalar> >::name()
01474       << "::takeSyncedStep("<<dt<<","<<toString(stepType)<<") ...\n"; 
01475   }
01476 
01477   return state_dt;
01478 
01479 }
01480 
01481 
01482 template<class Scalar> 
01483 Scalar ForwardSensitivityStepper<Scalar>::takeDecoupledStep(
01484   Scalar dt, StepSizeType stepType
01485   )
01486 {
01487 
01488 #ifdef ENABLE_RYTHMOS_TIMERS
01489   TEUCHOS_FUNC_TIME_MONITOR("Rythmos:ForwardSensitivityStepper::takeStep: decoupled");
01490 #endif
01491 
01492   using Teuchos::as;
01493   typedef Teuchos::ScalarTraits<Scalar> ST;
01494   typedef Teuchos::VerboseObjectTempState<InterpolationBufferBase<Scalar> > VOTSIBB;
01495   typedef Thyra::ModelEvaluatorBase MEB;
01496 
01497   RCP<Teuchos::FancyOStream> out = this->getOStream();
01498   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
01499   const bool lowTrace =
01500     ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) );
01501   const bool mediumTrace =
01502     ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_MEDIUM) );
01503   Teuchos::OSTab tab(out);
01504 
01505   if (lowTrace) {
01506     *out
01507       << "\nEntering " << TypeNameTraits<ForwardSensitivityStepper<Scalar> >::name()
01508       << "::takeDecoupledStep("<<dt<<","<<toString(stepType)<<") ...\n"; 
01509   }
01510   
01511   //
01512   // A) Take the sens timestep
01513   //
01514 
01515   if (lowTrace) {
01516     *out
01517       << "\nTaking sensitivity step using stepper : "
01518       << sensStepper_->description() << "\n";
01519   }
01520 
01521   Scalar sens_dt = -1.0;
01522   VOTSIBB sensStepper_outputTempState(sensStepper_,out,verbLevel);
01523   sens_dt = sensStepper_->takeStep(dt,stepType);
01524 
01525   if (mediumTrace) {
01526     const StepStatus<Scalar> sensStepStatus = sensStepper_->getStepStatus();
01527     *out << "\nSensitivity step status:\n" << sensStepStatus;
01528   }
01529 
01530   //
01531   // B) Wipe out all state interp buffer info before this sens timestep
01532   //
01533   
01534   //TEST_FOR_EXCEPT(true);
01535 
01536   if (lowTrace) {
01537     *out
01538       << "\nLeaving " << TypeNameTraits<ForwardSensitivityStepper<Scalar> >::name()
01539       << "::takeDecoupledStep("<<dt<<","<<toString(stepType)<<") ...\n"; 
01540   }
01541 
01542   return sens_dt;
01543   
01544 }
01545 
01546 
01547 } // namespace Rythmos
01548 
01549 
01550 #endif //RYTHMOS_FORWARD_SENSITIVITY_STEPPER_HPP
 All Classes Functions Variables Typedefs Friends