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_ForwardSensitivityModelEvaluatorBase.hpp"
00035 #include "Rythmos_ForwardSensitivityImplicitModelEvaluator.hpp"
00036 #include "Rythmos_ForwardSensitivityExplicitModelEvaluator.hpp"
00037 #include "Rythmos_StateAndForwardSensitivityModelEvaluator.hpp"
00038 #include "Rythmos_SolverAcceptingStepperBase.hpp"
00039 #include "Rythmos_IntegratorBase.hpp"
00040 #include "Rythmos_SingleResidualModelEvaluatorBase.hpp"
00041 #include "Thyra_ModelEvaluatorHelpers.hpp"
00042 #include "Thyra_LinearNonlinearSolver.hpp"
00043 #include "Thyra_ProductVectorBase.hpp"
00044 #include "Thyra_AssertOp.hpp"
00045 #include "Teuchos_ParameterListAcceptorDefaultBase.hpp"
00046 #include "Teuchos_VerboseObjectParameterListHelpers.hpp"
00047 #include "Teuchos_Assert.hpp"
00048 #include "Teuchos_as.hpp"
00049 
00050 
00051 namespace Rythmos {
00052 
00053 
00202 template<class Scalar> 
00203 class ForwardSensitivityStepper
00204   : virtual public StepperBase<Scalar>,
00205     virtual public Teuchos::ParameterListAcceptorDefaultBase
00206 {
00207 public:
00208 
00210   typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType ScalarMag;
00211 
00214 
00216   ForwardSensitivityStepper();
00217 
00271   void initializeSyncedSteppers(
00272     const RCP<const Thyra::ModelEvaluator<Scalar> > &stateModel,
00273     const int p_index,
00274     const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateBasePoint,
00275     const RCP<StepperBase<Scalar> > &stateStepper,
00276     const RCP<Thyra::NonlinearSolverBase<Scalar> > &stateTimeStepSolver,
00277     const RCP<StepperBase<Scalar> > &sensStepper = Teuchos::null,
00278     const RCP<Thyra::NonlinearSolverBase<Scalar> > &sensTimeStepSolver = Teuchos::null
00279     );
00280 
00316   void initializeDecoupledSteppers(
00317     const RCP<const Thyra::ModelEvaluator<Scalar> > &stateModel,
00318     const int p_index,
00319     const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateBasePoint,
00320     const RCP<StepperBase<Scalar> > &stateStepper,
00321     const RCP<Thyra::NonlinearSolverBase<Scalar> > &stateTimeStepSolver,
00322     const RCP<IntegratorBase<Scalar> > &stateIntegrator,
00323     const Scalar &finalTime,
00324     const RCP<StepperBase<Scalar> > &sensStepper = Teuchos::null,
00325     const RCP<Thyra::NonlinearSolverBase<Scalar> > &sensTimeStepSolver = Teuchos::null
00326     );
00327 
00331   RCP<const Thyra::ModelEvaluator<Scalar> >
00332   getStateModel() const;
00333 
00337   RCP<StepperBase<Scalar> >
00338   getNonconstStateStepper();
00339 
00343   RCP<const ForwardSensitivityModelEvaluatorBase<Scalar> >
00344   getFwdSensModel() const;
00345 
00352   RCP<const StateAndForwardSensitivityModelEvaluator<Scalar> >
00353   getStateAndFwdSensModel() const;
00354 
00356 
00359 
00361   void setParameterList(RCP<Teuchos::ParameterList> const& paramList);
00363   RCP<const Teuchos::ParameterList> getValidParameters() const;
00364 
00366 
00369 
00371   bool acceptsModel() const;
00372 
00374   void setModel(
00375     const RCP<const Thyra::ModelEvaluator<Scalar> > &model
00376     );
00377 
00384   RCP<const Thyra::ModelEvaluator<Scalar> >
00385   getModel() const;
00386 
00387   // RAB: ToDo: 2007/05/15: I need to talk with Todd about potentially
00388   // removing the setModel() and getModel() functions from the StepperBase
00389   // interface.  In the case of this forward sensitivity solver, I am not sure
00390   // that it makes a lot of sense to define a model.  This surely will not be
00391   // the model that a generic client would expect.  The assumption I am sure
00392   // would be that this model has the same space for x as the interpolation
00393   // buffer but that is not true in this case.
00394 
00408   void setInitialCondition(
00409     const Thyra::ModelEvaluatorBase::InArgs<Scalar> &state_and_sens_ic
00410     );
00411 
00413   Thyra::ModelEvaluatorBase::InArgs<Scalar> getInitialCondition() const;
00414 
00416   Scalar takeStep( Scalar dt, StepSizeType stepType );
00417 
00419   const StepStatus<Scalar> getStepStatus() const;
00420 
00422 
00425 
00432   RCP<const Thyra::VectorSpaceBase<Scalar> >
00433   get_x_space() const;
00434 
00436   void addPoints(
00437     const Array<Scalar>& time_vec,
00438     const Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >& x_vec,
00439     const Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >& xdot_vec
00440     );
00441 
00443   TimeRange<Scalar> getTimeRange() const;
00444 
00446   void getPoints(
00447     const Array<Scalar>& time_vec,
00448     Array<RCP<const Thyra::VectorBase<Scalar> > >* x_vec,
00449     Array<RCP<const Thyra::VectorBase<Scalar> > >* xdot_vec,
00450     Array<ScalarMag>* accuracy_vec
00451     ) const;
00452 
00454   void getNodes(Array<Scalar>* time_vec) const;
00455 
00457   void removeNodes(Array<Scalar>& time_vec);
00458 
00460   int getOrder() const;
00461 
00463 
00466 
00468   void initialize(
00469     const RCP<const Thyra::ModelEvaluator<Scalar> > &stateModel,
00470     const int p_index,
00471     const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateBasePoint,
00472     const RCP<StepperBase<Scalar> > &stateStepper,
00473     const RCP<Thyra::NonlinearSolverBase<Scalar> > &stateTimeStepSolver,
00474     const RCP<StepperBase<Scalar> > &sensStepper = Teuchos::null,
00475     const RCP<Thyra::NonlinearSolverBase<Scalar> > &sensTimeStepSolver = Teuchos::null
00476     )
00477     {
00478       initializeSyncedSteppers(
00479         stateModel, p_index, stateBasePoint, stateStepper, stateTimeStepSolver,
00480         sensStepper, sensTimeStepSolver
00481         );
00482     }
00483 
00485 
00486 private:
00487 
00488   // /////////////////////////
00489   // Private data members
00490 
00491   bool forceUpToDateW_;
00492   RCP<const Thyra::ModelEvaluator<Scalar> > stateModel_;
00493   Thyra::ModelEvaluatorBase::InArgs<Scalar> stateBasePoint_;
00494   RCP<StepperBase<Scalar> > stateStepper_;
00495   RCP<Thyra::NonlinearSolverBase<Scalar> > stateTimeStepSolver_;
00496   RCP<IntegratorBase<Scalar> > stateIntegrator_;
00497   Scalar finalTime_;
00498   Thyra::ModelEvaluatorBase::InArgs<Scalar> stateAndSensBasePoint_;
00499   RCP<StepperBase<Scalar> > sensStepper_;
00500   RCP<Thyra::NonlinearSolverBase<Scalar> > sensTimeStepSolver_;
00501 
00502   bool isSingleResidualStepper_;
00503   RCP<ForwardSensitivityModelEvaluatorBase<Scalar> > sensModel_;
00504   RCP<StateAndForwardSensitivityModelEvaluator<Scalar> > stateAndSensModel_;
00505   Thyra::ModelEvaluatorBase::InArgs<Scalar> stateBasePoint_t_;
00506 
00507   static const std::string forceUpToDateW_name_;
00508   static const bool forceUpToDateW_default_;
00509 
00510   // /////////////////////////
00511   // Private member functions
00512 
00513   void initializeCommon(
00514     const RCP<const Thyra::ModelEvaluator<Scalar> > &stateModel,
00515     const int p_index,
00516     const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateBasePoint,
00517     const RCP<StepperBase<Scalar> > &stateStepper,
00518     const RCP<Thyra::NonlinearSolverBase<Scalar> > &stateTimeStepSolver,
00519     const RCP<StepperBase<Scalar> > &sensStepper,
00520     const RCP<Thyra::NonlinearSolverBase<Scalar> > &sensTimeStepSolver
00521     );
00522 
00523   Scalar takeSyncedStep( Scalar dt, StepSizeType stepType );
00524 
00525   Scalar takeDecoupledStep( Scalar dt, StepSizeType stepType );
00526 
00527 };
00528 
00529 
00534 template<class Scalar> 
00535 inline
00536 RCP<ForwardSensitivityStepper<Scalar> >
00537 forwardSensitivityStepper()
00538 {
00539   return Teuchos::rcp(new ForwardSensitivityStepper<Scalar>());
00540 }
00541 
00542 
00547 template<class Scalar> 
00548 inline
00549 RCP<ForwardSensitivityStepper<Scalar> >
00550 forwardSensitivityStepper(
00551   const RCP<const Thyra::ModelEvaluator<Scalar> > &stateModel,
00552   const int p_index,
00553   const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateBasePoint,
00554   const RCP<StepperBase<Scalar> > &stateStepper,
00555   const RCP<Thyra::NonlinearSolverBase<Scalar> > &stateTimeStepSolver,
00556   const RCP<StepperBase<Scalar> > &sensStepper = Teuchos::null,
00557   const RCP<Thyra::NonlinearSolverBase<Scalar> > &sensTimeStepSolver = Teuchos::null
00558   )
00559 {
00560   RCP<ForwardSensitivityStepper<Scalar> >
00561     fwdSensStepper = Teuchos::rcp(new ForwardSensitivityStepper<Scalar>());
00562   fwdSensStepper->initialize(
00563     stateModel, p_index, stateBasePoint, stateStepper, stateTimeStepSolver );
00564   return fwdSensStepper;
00565 }
00566 
00567 
00573 template<class Scalar>
00574 int getParameterIndex(
00575   const ForwardSensitivityStepper<Scalar> &fwdSensStepper
00576   )
00577 {
00578   return fwdSensStepper.getFwdSensModel()->get_p_index();
00579 }
00580 
00581 
00588 template<class Scalar> 
00589 inline
00590 Thyra::ModelEvaluatorBase::InArgs<Scalar>
00591 createStateAndSensInitialCondition(
00592   const ForwardSensitivityStepper<Scalar> &fwdSensStepper,
00593   const Thyra::ModelEvaluatorBase::InArgs<Scalar> &state_ic
00594   )
00595 {
00596 
00597   using Teuchos::outArg;
00598   using Thyra::assign;
00599   typedef Thyra::ModelEvaluatorBase MEB;
00600   
00601   RCP<Thyra::VectorBase<Scalar> > s_bar_init
00602     = createMember(fwdSensStepper.getFwdSensModel()->get_x_space());
00603   assign( outArg(*s_bar_init), 0.0 );
00604   RCP<Thyra::VectorBase<Scalar> > s_bar_dot_init
00605     = createMember(fwdSensStepper.getFwdSensModel()->get_x_space());
00606   assign( outArg(*s_bar_dot_init), 0.0 );
00607   
00608   RCP<const Rythmos::StateAndForwardSensitivityModelEvaluator<Scalar> >
00609     stateAndSensModel = fwdSensStepper.getStateAndFwdSensModel();
00610   
00611   MEB::InArgs<Scalar>
00612     state_and_sens_ic = fwdSensStepper.getModel()->createInArgs();
00613   
00614   // Copy time, parameters etc.
00615   state_and_sens_ic.setArgs(state_ic);
00616   // Set initial condition for x_bar = [ x; s_bar ]
00617   state_and_sens_ic.set_x(
00618     stateAndSensModel->create_x_bar_vec(state_ic.get_x(), s_bar_init)
00619     );
00620   // Set initial condition for x_bar_dot = [ x_dot; s_bar_dot ]
00621   state_and_sens_ic.set_x_dot(
00622     stateAndSensModel->create_x_bar_vec(state_ic.get_x_dot(), s_bar_dot_init)
00623     );
00624 
00625   return state_and_sens_ic;
00626 
00627 }
00628 
00629 
00635 template<class Scalar> 
00636 inline
00637 Thyra::ModelEvaluatorBase::InArgs<Scalar>
00638 extractStateInitialCondition(
00639   const ForwardSensitivityStepper<Scalar> &fwdSensStepper,
00640   const Thyra::ModelEvaluatorBase::InArgs<Scalar> &state_and_sens_ic
00641   )
00642 {
00643 
00644   using Thyra::productVectorBase;
00645   typedef Thyra::ModelEvaluatorBase MEB;
00646 
00647   MEB::InArgs<Scalar>
00648     state_ic = fwdSensStepper.getStateModel()->createInArgs();
00649   
00650   // Copy time, parameters etc.
00651   state_ic.setArgs(state_and_sens_ic);
00652   state_ic.set_x(
00653     productVectorBase(state_and_sens_ic.get_x())->getVectorBlock(0));
00654   state_ic.set_x_dot(
00655     productVectorBase(state_and_sens_ic.get_x_dot())->getVectorBlock(0));
00656   
00657   return state_ic;
00658   
00659 }
00660 
00661 
00662 //
00663 // Implementation
00664 //
00665 
00666 
00667 // Static members
00668 
00669 
00670 template<class Scalar> 
00671 const std::string ForwardSensitivityStepper<Scalar>::forceUpToDateW_name_
00672 = "Force Up-To-Date Jacobian";
00673 
00674 template<class Scalar> 
00675 const bool ForwardSensitivityStepper<Scalar>::forceUpToDateW_default_
00676 = true;
00677 
00678 
00679 // Constructors, Intializers, Misc.
00680 
00681 
00682 template<class Scalar> 
00683 ForwardSensitivityStepper<Scalar>::ForwardSensitivityStepper()
00684   :forceUpToDateW_(false),
00685    isSingleResidualStepper_(false)
00686 {}
00687 
00688 
00689 template<class Scalar> 
00690 void ForwardSensitivityStepper<Scalar>::initializeSyncedSteppers(
00691   const RCP<const Thyra::ModelEvaluator<Scalar> > &stateModel,
00692   const int p_index,
00693   const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateBasePoint,
00694   const RCP<StepperBase<Scalar> > &stateStepper,
00695   const RCP<Thyra::NonlinearSolverBase<Scalar> > &stateTimeStepSolver,
00696   const RCP<StepperBase<Scalar> > &sensStepper,
00697   const RCP<Thyra::NonlinearSolverBase<Scalar> > &sensTimeStepSolver
00698   )
00699   
00700 {
00701   initializeCommon(stateModel, p_index, stateBasePoint, stateStepper,
00702     stateTimeStepSolver, sensStepper, sensTimeStepSolver );
00703 }
00704 
00705 
00706 template<class Scalar> 
00707 void ForwardSensitivityStepper<Scalar>::initializeDecoupledSteppers(
00708   const RCP<const Thyra::ModelEvaluator<Scalar> > &stateModel,
00709   const int p_index,
00710   const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateBasePoint,
00711   const RCP<StepperBase<Scalar> > &stateStepper,
00712   const RCP<Thyra::NonlinearSolverBase<Scalar> > &stateTimeStepSolver,
00713   const RCP<IntegratorBase<Scalar> > &stateIntegrator,
00714   const Scalar &finalTime,
00715   const RCP<StepperBase<Scalar> > &sensStepper,
00716   const RCP<Thyra::NonlinearSolverBase<Scalar> > &sensTimeStepSolver
00717   )
00718 {
00719   TEST_FOR_EXCEPT(is_null(stateIntegrator));
00720   initializeCommon(stateModel, p_index, stateBasePoint, stateStepper,
00721     stateTimeStepSolver, sensStepper, sensTimeStepSolver );
00722   stateIntegrator_ = stateIntegrator;
00723   finalTime_ = finalTime;
00724 }
00725 
00726   
00727 template<class Scalar> 
00728 RCP<const Thyra::ModelEvaluator<Scalar> >
00729 ForwardSensitivityStepper<Scalar>::getStateModel() const
00730 {
00731   return stateModel_;
00732 }
00733 
00734   
00735 template<class Scalar> 
00736 RCP<StepperBase<Scalar> >
00737 ForwardSensitivityStepper<Scalar>::getNonconstStateStepper()
00738 {
00739   return stateStepper_;
00740 }
00741 
00742   
00743 template<class Scalar> 
00744 RCP<const ForwardSensitivityModelEvaluatorBase<Scalar> >
00745 ForwardSensitivityStepper<Scalar>::getFwdSensModel() const
00746 {
00747   return sensModel_;
00748 }
00749 
00750   
00751 template<class Scalar> 
00752 RCP<const StateAndForwardSensitivityModelEvaluator<Scalar> >
00753 ForwardSensitivityStepper<Scalar>::getStateAndFwdSensModel() const
00754 {
00755   return stateAndSensModel_;
00756 }
00757 
00758 
00759 // Overridden from Teuchos::ParameterListAcceptor
00760 
00761 
00762 template<class Scalar> 
00763 void ForwardSensitivityStepper<Scalar>::setParameterList(
00764   RCP<Teuchos::ParameterList> const& paramList
00765   )
00766 {
00767   TEST_FOR_EXCEPT(is_null(paramList));
00768   paramList->validateParameters(*getValidParameters());
00769   this->setMyParamList(paramList);
00770   forceUpToDateW_ = paramList->get(forceUpToDateW_name_,forceUpToDateW_default_);
00771   Teuchos::readVerboseObjectSublist(&*paramList,this);
00772 }
00773 
00774 
00775 template<class Scalar> 
00776 RCP<const Teuchos::ParameterList>
00777 ForwardSensitivityStepper<Scalar>::getValidParameters() const
00778 {
00779   static RCP<const ParameterList> validPL;
00780   if (is_null(validPL) ) {
00781     RCP<ParameterList> pl = Teuchos::parameterList();
00782     pl->set( forceUpToDateW_name_, forceUpToDateW_default_,
00783       "If set to true, then the Jacobian matrix W used in the\n"
00784       "state timestep equation will be forced to be up to date\n"
00785       "with the final value for x for the nonlinear solve.  If\n"
00786       "you are willing to live with slightly less accurate sensitivities\n"
00787       "then set this to false."
00788       );
00789     Teuchos::setupVerboseObjectSublist(&*pl);
00790     validPL = pl;
00791   }
00792   return validPL;
00793 }
00794 
00795 
00796 // Overridden from StepperBase
00797 
00798 template<class Scalar> 
00799 bool ForwardSensitivityStepper<Scalar>::acceptsModel() const
00800 {
00801   return false;
00802 }
00803 
00804 template<class Scalar> 
00805 void ForwardSensitivityStepper<Scalar>::setModel(
00806   const RCP<const Thyra::ModelEvaluator<Scalar> > &model
00807   )
00808 {
00809   TEST_FOR_EXCEPT_MSG( true,
00810     "Error, this stepper subclass does not accept a model"
00811     " as defined by the StepperBase interface!");
00812 }
00813 
00814 
00815 template<class Scalar> 
00816 RCP<const Thyra::ModelEvaluator<Scalar> >
00817 ForwardSensitivityStepper<Scalar>::getModel() const
00818 {
00819   return stateAndSensModel_;
00820 }
00821 
00822 
00823 template<class Scalar> 
00824 void ForwardSensitivityStepper<Scalar>::setInitialCondition(
00825   const Thyra::ModelEvaluatorBase::InArgs<Scalar> &state_and_sens_ic
00826   )
00827 {
00828   
00829   typedef Thyra::ModelEvaluatorBase MEB;
00830 
00831   stateAndSensBasePoint_ = state_and_sens_ic;
00832 
00833   // Get the product vectors for x_bar = [ x; s_bar ] and x_bar_dot
00834 
00835   TEST_FOR_EXCEPTION(
00836     is_null(state_and_sens_ic.get_x()), std::logic_error,
00837     "Error, the initial condition for x_bar = [ x; s_bar ] can not be null!" );
00838 
00839   const RCP<const Thyra::ProductVectorBase<Scalar> >
00840     x_bar_init = Thyra::productVectorBase<Scalar>(
00841       state_and_sens_ic.get_x()
00842       );
00843 
00844   RCP<const Thyra::ProductVectorBase<Scalar> > x_bar_dot_init;
00845   if (state_and_sens_ic.supports(MEB::IN_ARG_x_dot)) {
00846       x_bar_dot_init = Thyra::productVectorBase<Scalar>(
00847         state_and_sens_ic.get_x_dot()
00848         );
00849   }
00850 
00851   // Remove x and x_dot from state_and_sens_ic_in to avoid cloning x and x dot!
00852   
00853   Thyra::ModelEvaluatorBase::InArgs<Scalar>
00854     state_and_sens_ic_no_x = state_and_sens_ic;
00855   state_and_sens_ic_no_x.set_x(Teuchos::null);
00856   if (state_and_sens_ic_no_x.supports(MEB::IN_ARG_x_dot)) {
00857     state_and_sens_ic_no_x.set_x_dot(Teuchos::null);
00858   }
00859 
00860   // Set initial condition for the state
00861 
00862   MEB::InArgs<Scalar> state_ic = stateModel_->createInArgs();
00863   state_ic.setArgs(state_and_sens_ic_no_x,true,true); // Set time, parameters etc.
00864   state_ic.set_x(x_bar_init->getVectorBlock(0)->clone_v());
00865   if (state_ic.supports(MEB::IN_ARG_x_dot)) {
00866     state_ic.set_x_dot(
00867         !is_null(x_bar_dot_init)
00868         ? x_bar_dot_init->getVectorBlock(0)->clone_v()
00869         : Teuchos::null
00870         );
00871   }
00872   stateStepper_->setInitialCondition(state_ic);
00873 
00874   // Set up the integrator if needed
00875   //if (!is_null(stateIntegrator_)) {
00876   //  stateIntegrator_->setStepper( stateStepper_, finalTime_ );
00877   //  sensModel_->setStateIntegrator( stateIntegrator_, state_ic );
00878   //}
00879 
00880   // Set initial condition for the sensitivities
00881   
00882   MEB::InArgs<Scalar> sens_ic = sensModel_->createInArgs();
00883   sens_ic.setArgs(state_and_sens_ic_no_x,true,true); // Set time etc.
00884   sens_ic.set_x(x_bar_init->getVectorBlock(1)->clone_v());
00885   if (sens_ic.supports(MEB::IN_ARG_x_dot)) {
00886     sens_ic.set_x_dot(
00887         !is_null(x_bar_dot_init)
00888         ? x_bar_dot_init->getVectorBlock(1)->clone_v()
00889         : Teuchos::null
00890         );
00891   }
00892   sensStepper_->setInitialCondition(sens_ic);
00893 
00894 }
00895 
00896 
00897 template<class Scalar> 
00898 Thyra::ModelEvaluatorBase::InArgs<Scalar>
00899 ForwardSensitivityStepper<Scalar>::getInitialCondition() const
00900 {
00901   return stateAndSensBasePoint_;
00902 }
00903  
00904 
00905 template<class Scalar> 
00906 Scalar
00907 ForwardSensitivityStepper<Scalar>::takeStep(
00908   Scalar dt, StepSizeType stepType
00909   )
00910 {
00911 
00912 #ifdef ENABLE_RYTHMOS_TIMERS
00913   TEUCHOS_FUNC_TIME_MONITOR("Rythmos:ForwardSensitivityStepper::takeStep");
00914 #endif
00915 
00916   if (!is_null(stateIntegrator_)) {
00917     return takeDecoupledStep(dt,stepType);
00918   }
00919 
00920   return takeSyncedStep(dt,stepType);
00921 
00922 }
00923 
00924 
00925 template<class Scalar> 
00926 const StepStatus<Scalar>
00927 ForwardSensitivityStepper<Scalar>::getStepStatus() const
00928 {
00929 
00930   const StepStatus<Scalar> sensStepStatus = sensStepper_->getStepStatus();
00931   StepStatus<Scalar> stepStatus;
00932   
00933   stepStatus.message = sensStepStatus.message;
00934   stepStatus.stepStatus = sensStepStatus.stepStatus;
00935   stepStatus.stepLETStatus = sensStepStatus.stepLETStatus;
00936   stepStatus.stepSize = sensStepStatus.stepSize;
00937   stepStatus.order = sensStepStatus.order;
00938   stepStatus.time = sensStepStatus.time;
00939   stepStatus.stepLETValue = sensStepStatus.stepLETValue;
00940   stepStatus.extraParameters = sensStepStatus.extraParameters;
00941 
00942   if (is_null(stateIntegrator_)) {
00943     const StepStatus<Scalar> 
00944       stateStepStatus = stateStepper_->getStepStatus();
00945     if (!is_null(stateStepStatus.solution) && !is_null(sensStepStatus.solution))
00946       stepStatus.solution = stateAndSensModel_->create_x_bar_vec(
00947         stateStepStatus.solution, sensStepStatus.solution
00948         );
00949     if (!is_null(stateStepStatus.solutionDot) && !is_null(sensStepStatus.solutionDot))
00950       stepStatus.solutionDot = stateAndSensModel_->create_x_bar_vec(
00951         stateStepStatus.solutionDot, sensStepStatus.solutionDot
00952         );
00953   }
00954 
00955   return stepStatus;
00956 
00957 }
00958 
00959 
00960 // Overridden from InterpolationBufferBase
00961 
00962 
00963 template<class Scalar> 
00964 RCP<const Thyra::VectorSpaceBase<Scalar> >
00965 ForwardSensitivityStepper<Scalar>::get_x_space() const
00966 {
00967   return stateAndSensModel_->get_x_space();
00968 }
00969 
00970 
00971 template<class Scalar> 
00972 void ForwardSensitivityStepper<Scalar>::addPoints(
00973   const Array<Scalar>& time_vec,
00974   const Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >& x_vec,
00975   const Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >& xdot_vec
00976   )
00977 {
00978   TEST_FOR_EXCEPT("Not implemented addPoints(...) yet but we could if we wanted!");
00979 }
00980 
00981 
00982 template<class Scalar> 
00983 TimeRange<Scalar>
00984 ForwardSensitivityStepper<Scalar>::getTimeRange() const
00985 {
00986   return sensStepper_->getTimeRange();
00987 }
00988 
00989 
00990 template<class Scalar> 
00991 void ForwardSensitivityStepper<Scalar>::getPoints(
00992   const Array<Scalar>& time_vec,
00993   Array<RCP<const Thyra::VectorBase<Scalar> > >* x_bar_vec,
00994   Array<RCP<const Thyra::VectorBase<Scalar> > >* x_bar_dot_vec,
00995   Array<ScalarMag>* accuracy_vec
00996   ) const
00997 {
00998 
00999   using Teuchos::as;
01000 
01001 #ifdef RYTHMOS_DEBUG
01002   TEST_FOR_EXCEPT( as<int>(time_vec.size()) == 0 );
01003 #endif
01004 
01005   const int numTimePoints = time_vec.size();
01006 
01007   if (x_bar_vec)
01008     x_bar_vec->clear();
01009 
01010   if (x_bar_dot_vec)
01011     x_bar_dot_vec->clear();
01012   
01013   Array<RCP<const Thyra::VectorBase<Scalar> > >
01014     x_vec, x_dot_vec;
01015 
01016   if (!is_null(stateIntegrator_)) {
01017     stateIntegrator_->getPoints(
01018       time_vec,
01019       x_bar_vec ? &x_vec: 0,
01020       x_bar_dot_vec ? &x_dot_vec: 0,
01021       0 // Ignoring accuracy from state for now!
01022       );
01023   }
01024   else {
01025     stateStepper_->getPoints(
01026       time_vec,
01027       x_bar_vec ? &x_vec: 0,
01028       x_bar_dot_vec ? &x_dot_vec: 0,
01029       0 // Ignoring accuracy from state for now!
01030       );
01031   }
01032   
01033   Array<RCP<const Thyra::VectorBase<Scalar> > >
01034     s_bar_vec, s_bar_dot_vec;
01035 
01036   sensStepper_->getPoints(
01037     time_vec,
01038     x_bar_vec ? &s_bar_vec: 0,
01039     x_bar_dot_vec ? &s_bar_dot_vec: 0,
01040     accuracy_vec
01041     );
01042   
01043   if ( x_bar_vec ) {
01044     for ( int i = 0; i < numTimePoints; ++i ) {
01045       x_bar_vec->push_back(
01046         stateAndSensModel_->create_x_bar_vec(x_vec[i],s_bar_vec[i])
01047         );
01048     }
01049   }
01050   
01051   if ( x_bar_dot_vec ) {
01052     for ( int i = 0; i < numTimePoints; ++i ) {
01053       x_bar_dot_vec->push_back(
01054         stateAndSensModel_->create_x_bar_vec(x_dot_vec[i],s_bar_dot_vec[i])
01055         );
01056     }
01057   }
01058 
01059 }
01060 
01061 
01062 template<class Scalar>
01063 void ForwardSensitivityStepper<Scalar>::getNodes(
01064   Array<Scalar>* time_vec
01065   ) const
01066 {
01067   TEUCHOS_ASSERT( time_vec != NULL );
01068   time_vec->clear();
01069   if (is_null(stateIntegrator_) && is_null(stateStepper_)) {
01070     return;
01071   }
01072   if (!is_null(stateIntegrator_)) {
01073     stateIntegrator_->getNodes(time_vec);
01074   }
01075   else {
01076     stateStepper_->getNodes(time_vec);
01077   }
01078 }
01079 
01080 
01081 template<class Scalar> 
01082 void ForwardSensitivityStepper<Scalar>::removeNodes(
01083   Array<Scalar>& time_vec
01084   )
01085 {
01086   TEST_FOR_EXCEPT("Not implemented yet but we can!");
01087 }
01088 
01089 
01090 template<class Scalar> 
01091 int ForwardSensitivityStepper<Scalar>::getOrder() const
01092 {
01093   return sensStepper_->getOrder();
01094   // Note: This assumes that stateStepper will have the same order!
01095 }
01096 
01097 
01098 // private
01099 
01100 
01101 template<class Scalar> 
01102 void ForwardSensitivityStepper<Scalar>::initializeCommon(
01103   const RCP<const Thyra::ModelEvaluator<Scalar> > &stateModel,
01104   const int p_index,
01105   const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateBasePoint,
01106   const RCP<StepperBase<Scalar> > &stateStepper,
01107   const RCP<Thyra::NonlinearSolverBase<Scalar> > &stateTimeStepSolver,
01108   const RCP<StepperBase<Scalar> > &sensStepper,
01109   const RCP<Thyra::NonlinearSolverBase<Scalar> > &sensTimeStepSolver
01110   )
01111 {
01112 
01113   using Teuchos::rcp_dynamic_cast;
01114 
01115   typedef Thyra::ModelEvaluatorBase MEB;
01116 
01117   //
01118   // Validate input
01119   //
01120 
01121   TEST_FOR_EXCEPT( is_null(stateModel) );
01122   TEST_FOR_EXCEPT( is_null(stateStepper) );
01123   if (stateStepper->isImplicit()) {
01124     TEST_FOR_EXCEPT( is_null(stateTimeStepSolver) ); // allow to be null for explicit methods
01125   }
01126 
01127 
01128   //
01129   // Create the sensModel which will do some more validation
01130   //
01131   
01132   RCP<ForwardSensitivityModelEvaluatorBase<Scalar> > sensModel;
01133   MEB::InArgs<Scalar> stateModelInArgs = stateModel->createInArgs();
01134   if (stateModelInArgs.supports(MEB::IN_ARG_x_dot)) {
01135     // Implicit DE formulation
01136     sensModel = Teuchos::rcp(new ForwardSensitivityImplicitModelEvaluator<Scalar>);
01137   }
01138   else {
01139     // Explicit DE formulation
01140     sensModel = Teuchos::rcp(new ForwardSensitivityExplicitModelEvaluator<Scalar>);
01141   }
01142 
01143   sensModel->initializeStructure(stateModel,p_index);
01144   
01145   //
01146   // Get the input objects
01147   //
01148 
01149   stateModel_ = stateModel;
01150 
01151   stateBasePoint_ = stateBasePoint;
01152 
01153   stateStepper_ = stateStepper;
01154   
01155   stateTimeStepSolver_ = stateTimeStepSolver;
01156 
01157   sensModel_ = sensModel;
01158 
01159   stateAndSensModel_ = Teuchos::rcp(new StateAndForwardSensitivityModelEvaluator<Scalar>);
01160   stateAndSensModel_->initializeStructure(sensModel_);
01161 
01162   if (!is_null(sensStepper)) {
01163     sensStepper_ = sensStepper;
01164   }
01165   else {
01166     sensStepper_ = stateStepper_->cloneStepperAlgorithm();
01167     TEST_FOR_EXCEPTION(
01168       is_null(sensStepper_), std::logic_error,
01169       "Error, if the client does not pass in a stepper for the senitivity\n"
01170       "equations then the stateStepper object must support cloning to create\n"
01171       "the sensitivity stepper!"
01172       );
01173   }
01174 
01175   if (!is_null(sensTimeStepSolver)) {
01176     sensTimeStepSolver_ = sensTimeStepSolver;
01177   }
01178   else {
01179     RCP<Thyra::LinearNonlinearSolver<Scalar> >
01180       linearNonlinearSolver(new Thyra::LinearNonlinearSolver<Scalar>);
01181     // ToDo: Set tolerance on the nonlinear solver???
01182     sensTimeStepSolver_ = linearNonlinearSolver;
01183   }
01184 
01185   //
01186   // Setup the steppers
01187   //
01188 
01189   isSingleResidualStepper_ = true; // ToDo: Add dynamic cast on
01190                                    // stateTimeStepSolver to check this!
01191 
01192   stateStepper_->setModel(stateModel_);
01193   if (stateStepper_->isImplicit()) {
01194     rcp_dynamic_cast<SolverAcceptingStepperBase<Scalar> >(
01195         stateStepper_,true)->setSolver(stateTimeStepSolver_);
01196   }
01197   sensStepper_->setModel(sensModel_);
01198   if (sensStepper_->isImplicit()) {
01199     rcp_dynamic_cast<SolverAcceptingStepperBase<Scalar> >(
01200         sensStepper_,true)->setSolver(sensTimeStepSolver_);
01201   }
01202 
01203   stateBasePoint_t_ = stateModel_->createInArgs();
01204 
01205   // 2007/05/18: rabartl: ToDo: Move the above initialization code to give
01206   // setInitializeCondition(...) a chance to set the initial condition.
01207 
01208 }
01209 
01210 
01211 template<class Scalar> 
01212 Scalar ForwardSensitivityStepper<Scalar>::takeSyncedStep(
01213   Scalar dt, StepSizeType stepType
01214   )
01215 {
01216 
01217 #ifdef ENABLE_RYTHMOS_TIMERS
01218   TEUCHOS_FUNC_TIME_MONITOR("Rythmos:ForwardSensitivityStepper::takeStep: synced");
01219 #endif
01220 
01221   using Teuchos::as;
01222   typedef Teuchos::ScalarTraits<Scalar> ST;
01223   typedef Teuchos::VerboseObjectTempState<InterpolationBufferBase<Scalar> > VOTSIBB;
01224   typedef Thyra::ModelEvaluatorBase MEB;
01225 
01226   RCP<Teuchos::FancyOStream> out = this->getOStream();
01227   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
01228   const bool lowTrace =
01229     ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) );
01230   const bool mediumTrace =
01231     ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_MEDIUM) );
01232   Teuchos::OSTab tab(out);
01233 
01234   if (lowTrace) {
01235     *out
01236       << "\nEntering " << TypeNameTraits<ForwardSensitivityStepper<Scalar> >::name()
01237       << "::takeSyncedStep("<<dt<<","<<toString(stepType)<<") ...\n"; 
01238   }
01239 
01240   //
01241   // A) Compute the state timestep
01242   //
01243 
01244   if (lowTrace) {
01245     *out
01246       << "\nTaking state step using stepper : "
01247       << stateStepper_->description() << "\n";
01248   }
01249 
01250   Scalar state_dt = -1.0;
01251   {
01252 #ifdef ENABLE_RYTHMOS_TIMERS
01253     TEUCHOS_FUNC_TIME_MONITOR("Rythmos:ForwardSensitivityStepper::takeStep: stateStep");
01254 #endif
01255     VOTSIBB stateStepper_outputTempState(stateStepper_,out,verbLevel);
01256     state_dt = stateStepper_->takeStep(dt,stepType);
01257   }
01258 
01259   if (state_dt < Scalar(-ST::one())) {
01260     if (lowTrace)
01261       *out << "\nThe state stepper has failed so return a failed timestep!\n";
01262     return state_dt;
01263   }
01264   
01265   const StepStatus<Scalar> stateStepStatus = stateStepper_->getStepStatus();
01266   if (mediumTrace)
01267     *out << "\nState step status:\n" << stateStepStatus;
01268 
01269   //
01270   // B) Setup the sensitivity model for this timestep
01271   //
01272 
01273   {
01274 
01275 #ifdef ENABLE_RYTHMOS_TIMERS
01276     TEUCHOS_FUNC_TIME_MONITOR("Rythmos:ForwardSensitivityStepper::takeStep: updateSensModel");
01277 #endif
01278     
01279     TEST_FOR_EXCEPTION(
01280       stateStepStatus.stepStatus != STEP_STATUS_CONVERGED, std::logic_error,
01281       "Error, the status should be converged since a positive step size was returned!"
01282       );
01283     
01284     const Scalar curr_t = stateStepStatus.time;
01285 
01286     RCP<ForwardSensitivityImplicitModelEvaluator<Scalar> > implicitSensModel = 
01287       Teuchos::rcp_dynamic_cast<ForwardSensitivityImplicitModelEvaluator<Scalar> >(sensModel_,false);
01288     if (!is_null(implicitSensModel)) { // ForwardSensitivityImplicitModelEvaluator
01289       // Get both x and x_dot since these are needed compute other derivative
01290       // objects at these points.
01291       RCP<const Thyra::VectorBase<Scalar> > x, x_dot;
01292       get_x_and_x_dot(*stateStepper_,curr_t,&x,&x_dot);
01293       
01294       stateBasePoint_t_ = stateBasePoint_;
01295       stateBasePoint_t_.set_x_dot( x_dot );
01296       stateBasePoint_t_.set_x( x );
01297       stateBasePoint_t_.set_t( curr_t );
01298 
01299       // Grab the SingleResidualModel that was used to compute the state timestep.
01300       // From this, we can get the constants that where used to compute W!
01301       RCP<const Rythmos::SingleResidualModelEvaluatorBase<Scalar> >
01302         singleResidualModel
01303         = Teuchos::rcp_dynamic_cast<const Rythmos::SingleResidualModelEvaluatorBase<Scalar> >(
01304           stateTimeStepSolver_->getModel(), true
01305           );
01306       const Scalar
01307         coeff_x_dot = singleResidualModel->get_coeff_x_dot(),
01308         coeff_x = singleResidualModel->get_coeff_x();
01309       
01310       // Get W (and force an update if not up to date already)
01311 
01312       if (mediumTrace && forceUpToDateW_)
01313         *out << "\nForcing an update of W at the converged state timestep ...\n";
01314       
01315       RCP<Thyra::LinearOpWithSolveBase<Scalar> >
01316         W_tilde = stateTimeStepSolver_->get_nonconst_W(forceUpToDateW_);
01317       
01318       TEST_FOR_EXCEPTION(
01319         is_null(W_tilde), std::logic_error,
01320         "Error, the W from the state time step be non-null!"
01321         );
01322       
01323       implicitSensModel->initializeState( stateBasePoint_t_, W_tilde, coeff_x_dot, coeff_x );
01324 
01325     } else { // explicit ME:
01326       // ForwardSensitivityExplicitModelEvaluator
01327       RCP<ForwardSensitivityExplicitModelEvaluator<Scalar> > explicitSensModel = 
01328         Teuchos::rcp_dynamic_cast<ForwardSensitivityExplicitModelEvaluator<Scalar> > (sensModel_,true);
01329       
01330       RCP<const Thyra::VectorBase<Scalar> > x;
01331       x = get_x(*stateStepper_,curr_t);
01332       
01333       stateBasePoint_t_ = stateBasePoint_;
01334       stateBasePoint_t_.set_x( x );
01335       stateBasePoint_t_.set_t( curr_t );
01336 
01337       explicitSensModel->initializePointState( stateBasePoint_t_ );
01338     }
01339 
01340   }
01341   // 2007/09/04: rabartl: ToDo: Move the above code that initializes the
01342   // senstivity stepper from the state stepper as a strategy object that can
01343   // be redefined.  This will allow me to make this class 100% general for all
01344   // time stepper methods ... Including explicit steppers!
01345 
01346   //
01347   // C) Compute the sensitivity timestep for the exact same timestep as was
01348   // used for the state solve.
01349   //
01350 
01351   if (lowTrace) {
01352     *out
01353       << "\nTaking sensitivity step using stepper : "
01354       << sensStepper_->description() << "\n";
01355   }
01356 
01357   Scalar sens_dt = -1.0;
01358   {
01359 #ifdef ENABLE_RYTHMOS_TIMERS
01360     TEUCHOS_FUNC_TIME_MONITOR("Rythmos:ForwardSensitivityStepper::takeStep: sensStep");
01361 #endif
01362     // Copy the step control data to make sure that the sensStepper takes the
01363     // same type of step that the statStepper took.  This is needed to ensure
01364     // that the W matrix is the same for one.
01365     sensStepper_->setStepControlData(*stateStepper_);
01366     VOTSIBB sensStepper_outputTempState(sensStepper_,out,verbLevel);
01367     sens_dt = sensStepper_->takeStep(state_dt,STEP_TYPE_FIXED);
01368   }
01369 
01370   if (mediumTrace) {
01371     const StepStatus<Scalar> sensStepStatus = sensStepper_->getStepStatus();
01372     *out << "\nSensitivity step status:\n" << sensStepStatus;
01373   }
01374   
01375   TEST_FOR_EXCEPTION(
01376     sens_dt != state_dt, std::logic_error,
01377     "Error, the sensitivity step failed for some reason.  We should\n"
01378     "just return a negative step size and reject the step but currently\n"
01379     "there is no way to roll back the state timestep it for back to\n"
01380     "the status before this function was called!"
01381     );
01382 
01383   // 2007/05/18: rabartl: ToDo: If stepType == STEP_TYPE_VARIABLE and the state
01384   // timestep sucessed but the sensitivity timestep failed, then we need to
01385   // really throw an excpetion because there is nothing that we can really do
01386   // here!
01387 
01388   // 2007/05/18: rabartl: ToDo: Replace the above std::logic_error type with
01389   // a Rythmos::CatastrophicFailure or just use Thyra::CatastrophicFailure!
01390 
01391   if (lowTrace) {
01392     *out
01393       << "\nLeaving " << TypeNameTraits<ForwardSensitivityStepper<Scalar> >::name()
01394       << "::takeSyncedStep("<<dt<<","<<toString(stepType)<<") ...\n"; 
01395   }
01396 
01397   return state_dt;
01398 
01399 }
01400 
01401 
01402 template<class Scalar> 
01403 Scalar ForwardSensitivityStepper<Scalar>::takeDecoupledStep(
01404   Scalar dt, StepSizeType stepType
01405   )
01406 {
01407 
01408 #ifdef ENABLE_RYTHMOS_TIMERS
01409   TEUCHOS_FUNC_TIME_MONITOR("Rythmos:ForwardSensitivityStepper::takeStep: decoupled");
01410 #endif
01411 
01412   using Teuchos::as;
01413   typedef Teuchos::ScalarTraits<Scalar> ST;
01414   typedef Teuchos::VerboseObjectTempState<InterpolationBufferBase<Scalar> > VOTSIBB;
01415   typedef Thyra::ModelEvaluatorBase MEB;
01416 
01417   RCP<Teuchos::FancyOStream> out = this->getOStream();
01418   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
01419   const bool lowTrace =
01420     ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) );
01421   const bool mediumTrace =
01422     ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_MEDIUM) );
01423   Teuchos::OSTab tab(out);
01424 
01425   if (lowTrace) {
01426     *out
01427       << "\nEntering " << TypeNameTraits<ForwardSensitivityStepper<Scalar> >::name()
01428       << "::takeDecoupledStep("<<dt<<","<<toString(stepType)<<") ...\n"; 
01429   }
01430   
01431   //
01432   // A) Take the sens timestep
01433   //
01434 
01435   if (lowTrace) {
01436     *out
01437       << "\nTaking sensitivity step using stepper : "
01438       << sensStepper_->description() << "\n";
01439   }
01440 
01441   Scalar sens_dt = -1.0;
01442   VOTSIBB sensStepper_outputTempState(sensStepper_,out,verbLevel);
01443   sens_dt = sensStepper_->takeStep(dt,stepType);
01444 
01445   if (mediumTrace) {
01446     const StepStatus<Scalar> sensStepStatus = sensStepper_->getStepStatus();
01447     *out << "\nSensitivity step status:\n" << sensStepStatus;
01448   }
01449 
01450   //
01451   // B) Wipe out all state interp buffer info before this sens timestep
01452   //
01453   
01454   //TEST_FOR_EXCEPT(true);
01455 
01456   if (lowTrace) {
01457     *out
01458       << "\nLeaving " << TypeNameTraits<ForwardSensitivityStepper<Scalar> >::name()
01459       << "::takeDecoupledStep("<<dt<<","<<toString(stepType)<<") ...\n"; 
01460   }
01461 
01462   return sens_dt;
01463   
01464 }
01465 
01466 
01467 } // namespace Rythmos
01468 
01469 
01470 #endif //RYTHMOS_FORWARD_SENSITIVITY_STEPPER_HPP

Generated on Wed May 12 21:25:43 2010 for Rythmos - Transient Integration for Differential Equations by  doxygen 1.4.7