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_ForwardSensitivityModelEvaluator.hpp"
00035 #include "Rythmos_StateAndForwardSensitivityModelEvaluator.hpp"
00036 #include "Rythmos_SolverAcceptingStepperBase.hpp"
00037 #include "Rythmos_SingleResidualModelEvaluatorBase.hpp"
00038 #include "Thyra_ModelEvaluatorHelpers.hpp"
00039 #include "Thyra_LinearNonlinearSolver.hpp"
00040 #include "Thyra_AssertOp.hpp"
00041 #include "Teuchos_ParameterListAcceptorDefaultBase.hpp"
00042 #include "Teuchos_VerboseObjectParameterListHelpers.hpp"
00043 #include "Teuchos_Assert.hpp"
00044 #include "Teuchos_as.hpp"
00045 
00046 
00047 namespace Rythmos {
00048 
00049 
00198 template<class Scalar> 
00199 class ForwardSensitivityStepper
00200   : virtual public StepperBase<Scalar>,
00201     virtual public Teuchos::ParameterListAcceptorDefaultBase
00202 {
00203 public:
00204 
00206   typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType ScalarMag;
00207 
00210 
00212   ForwardSensitivityStepper();
00213 
00266   void initialize(
00267     const RCP<const Thyra::ModelEvaluator<Scalar> > &stateModel,
00268     const int p_index,
00269     const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateBasePoint,
00270     const RCP<StepperBase<Scalar> > &stateStepper,
00271     const RCP<Thyra::NonlinearSolverBase<Scalar> > &stateTimeStepSolver,
00272     const RCP<StepperBase<Scalar> > &sensStepper = Teuchos::null,
00273     const RCP<Thyra::NonlinearSolverBase<Scalar> > &sensTimeStepSolver = Teuchos::null
00274     );
00275 
00279   RCP<const Thyra::ModelEvaluator<Scalar> >
00280   getStateModel() const;
00281 
00285   RCP<const ForwardSensitivityModelEvaluator<Scalar> >
00286   getFwdSensModel() const;
00287 
00294   RCP<const StateAndForwardSensitivityModelEvaluator<Scalar> >
00295   getStateAndFwdSensModel() const;
00296 
00298 
00301 
00303   void setParameterList(RCP<Teuchos::ParameterList> const& paramList);
00305   RCP<const Teuchos::ParameterList> getValidParameters() const;
00306 
00308 
00311 
00313   bool acceptsModel() const;
00314 
00316   void setModel(
00317     const RCP<const Thyra::ModelEvaluator<Scalar> > &model
00318     );
00319 
00326   RCP<const Thyra::ModelEvaluator<Scalar> >
00327   getModel() const;
00328 
00329   // RAB: ToDo: 2007/05/15: I need to talk with Todd about potentially
00330   // removing the setModel() and getModel() functions from the StepperBase
00331   // interface.  In the case of this forward sensitivity solver, I am not sure
00332   // that it makes a lot of sense to define a model.  This surely will not be
00333   // the model that a generic client would expect.  The assumption I am sure
00334   // would be that this model has the same space for x as the interpolation
00335   // buffer but that is not true in this case.
00336 
00350   void setInitialCondition(
00351     const Thyra::ModelEvaluatorBase::InArgs<Scalar> &state_and_sens_ic
00352     );
00353 
00355   Scalar takeStep( Scalar dt, StepSizeType stepType );
00356 
00358   const StepStatus<Scalar> getStepStatus() const;
00359 
00361 
00364 
00371   RCP<const Thyra::VectorSpaceBase<Scalar> >
00372   get_x_space() const;
00373 
00375   void addPoints(
00376     const Array<Scalar>& time_vec,
00377     const Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >& x_vec,
00378     const Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >& xdot_vec
00379     );
00380 
00382   TimeRange<Scalar> getTimeRange() const;
00383 
00385   void getPoints(
00386     const Array<Scalar>& time_vec,
00387     Array<RCP<const Thyra::VectorBase<Scalar> > >* x_vec,
00388     Array<RCP<const Thyra::VectorBase<Scalar> > >* xdot_vec,
00389     Array<ScalarMag>* accuracy_vec
00390     ) const;
00391 
00393   void getNodes(Array<Scalar>* time_vec) const;
00394 
00396   void removeNodes(Array<Scalar>& time_vec);
00397 
00399   int getOrder() const;
00400 
00402 
00403 private:
00404 
00405   // /////////////////////////
00406   // Private data members
00407 
00408   bool forceUpToDateW_;
00409   RCP<const Thyra::ModelEvaluator<Scalar> > stateModel_;
00410   Thyra::ModelEvaluatorBase::InArgs<Scalar> stateBasePoint_;
00411   RCP<StepperBase<Scalar> > stateStepper_;
00412   RCP<Thyra::NonlinearSolverBase<Scalar> > stateTimeStepSolver_;
00413   RCP<StepperBase<Scalar> > sensStepper_;
00414   RCP<Thyra::NonlinearSolverBase<Scalar> > sensTimeStepSolver_;
00415 
00416   bool isSingleResidualStepper_;
00417   RCP<ForwardSensitivityModelEvaluator<Scalar> > sensModel_;
00418   RCP<StateAndForwardSensitivityModelEvaluator<Scalar> > stateAndSensModel_;
00419   Thyra::ModelEvaluatorBase::InArgs<Scalar> stateBasePoint_t_;
00420 
00421   static const std::string forceUpToDateW_name_;
00422   static const bool forceUpToDateW_default_;
00423 
00424 };
00425 
00426 
00431 template<class Scalar> 
00432 inline
00433 RCP<ForwardSensitivityStepper<Scalar> >
00434 forwardSensitivityStepper(
00435   const RCP<const Thyra::ModelEvaluator<Scalar> > &stateModel,
00436   const int p_index,
00437   const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateBasePoint,
00438   const RCP<StepperBase<Scalar> > &stateStepper,
00439   const RCP<Thyra::NonlinearSolverBase<Scalar> > &stateTimeStepSolver,
00440   const RCP<StepperBase<Scalar> > &sensStepper = Teuchos::null,
00441   const RCP<Thyra::NonlinearSolverBase<Scalar> > &sensTimeStepSolver = Teuchos::null
00442   )
00443 {
00444   RCP<ForwardSensitivityStepper<Scalar> >
00445     fwdSensStepper = Teuchos::rcp(new ForwardSensitivityStepper<Scalar>());
00446   fwdSensStepper->initialize(
00447     stateModel, p_index, stateBasePoint, stateStepper, stateTimeStepSolver );
00448   return fwdSensStepper;
00449 }
00450 
00451 
00457 template<class Scalar>
00458 int getParameterIndex(
00459   const ForwardSensitivityStepper<Scalar> &fwdSensStepper
00460   )
00461 {
00462   return fwdSensStepper.getFwdSensModel()->get_p_index();
00463 }
00464 
00465 
00466 //
00467 // Implementation
00468 //
00469 
00470 
00471 // Static members
00472 
00473 
00474 template<class Scalar> 
00475 const std::string ForwardSensitivityStepper<Scalar>::forceUpToDateW_name_
00476 = "Force Up-To-Date Jacobian";
00477 
00478 template<class Scalar> 
00479 const bool ForwardSensitivityStepper<Scalar>::forceUpToDateW_default_
00480 = true;
00481 
00482 // Constructors, Intializers, Misc.
00483 
00484 
00485 template<class Scalar> 
00486 ForwardSensitivityStepper<Scalar>::ForwardSensitivityStepper()
00487   :forceUpToDateW_(false),
00488    isSingleResidualStepper_(false)
00489 {}
00490 
00491 
00492 template<class Scalar> 
00493 void ForwardSensitivityStepper<Scalar>::initialize(
00494   const RCP<const Thyra::ModelEvaluator<Scalar> > &stateModel,
00495   const int p_index,
00496   const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateBasePoint,
00497   const RCP<StepperBase<Scalar> > &stateStepper,
00498   const RCP<Thyra::NonlinearSolverBase<Scalar> > &stateTimeStepSolver,
00499   const RCP<StepperBase<Scalar> > &sensStepper,
00500   const RCP<Thyra::NonlinearSolverBase<Scalar> > &sensTimeStepSolver
00501   )
00502   
00503 {
00504 
00505   using Teuchos::rcp_dynamic_cast;
00506   using Teuchos::tuple;
00507 
00508   typedef Thyra::ModelEvaluatorBase MEB;
00509 
00510   //
00511   // Validate input
00512   //
00513 
00514   TEST_FOR_EXCEPT( is_null(stateModel) );
00515   TEST_FOR_EXCEPT( is_null(stateStepper) );
00516   TEST_FOR_EXCEPT( is_null(stateTimeStepSolver) ); // ToDo: allow to be null for explicit methods!
00517 
00518 
00519   //
00520   // Create the sensModel which will do some more validation
00521   //
00522   
00523   RCP<ForwardSensitivityModelEvaluator<Scalar> >
00524     sensModel = Teuchos::rcp(new ForwardSensitivityModelEvaluator<Scalar>);
00525   sensModel->initializeStructure(stateModel,p_index);
00526   
00527   //
00528   // Get the input objects
00529   //
00530 
00531   stateModel_ = stateModel;
00532 
00533   stateBasePoint_ = stateBasePoint;
00534 
00535   stateStepper_ = stateStepper;
00536   
00537   stateTimeStepSolver_ = stateTimeStepSolver;
00538 
00539   sensModel_ = sensModel;
00540 
00541   stateAndSensModel_ = Teuchos::rcp(new StateAndForwardSensitivityModelEvaluator<Scalar>);
00542   stateAndSensModel_->initializeStructure(sensModel_);
00543 
00544   if (!is_null(sensStepper)) {
00545     sensStepper_ = sensStepper;
00546   }
00547   else {
00548     sensStepper_ = stateStepper_->cloneStepperAlgorithm();
00549     TEST_FOR_EXCEPTION(
00550       is_null(sensStepper_), std::logic_error,
00551       "Error, if the client does not pass in a stepper for the senitivity\n"
00552       "equations then the stateStepper object must support cloning to create\n"
00553       "the sensitivity stepper!"
00554       );
00555   }
00556 
00557   if (!is_null(sensTimeStepSolver)) {
00558     sensTimeStepSolver_ = sensTimeStepSolver;
00559   }
00560   else {
00561     RCP<Thyra::LinearNonlinearSolver<Scalar> >
00562       linearNonlinearSolver(new Thyra::LinearNonlinearSolver<Scalar>);
00563     // ToDo: Set tolerance on the nonlinear solver???
00564     sensTimeStepSolver_ = linearNonlinearSolver;
00565   }
00566 
00567   //
00568   // Setup the steppers
00569   //
00570 
00571   isSingleResidualStepper_ = true; // ToDo: Add dynamic cast on
00572                                    // stateTimeStepSolver to check this!
00573 
00574   stateStepper_->setModel(stateModel_);
00575   rcp_dynamic_cast<SolverAcceptingStepperBase<Scalar> >(
00576     stateStepper_,true)->setSolver(stateTimeStepSolver_);
00577   sensStepper_->setModel(sensModel_);
00578   rcp_dynamic_cast<SolverAcceptingStepperBase<Scalar> >(
00579     sensStepper_,true)->setSolver(sensTimeStepSolver_);
00580 
00581   stateBasePoint_t_ = stateModel_->createInArgs();
00582 
00583   // 2007/05/18: rabartl: ToDo: Move the above initialization code to give
00584   // setInitializeCondition(...) a chance to set the initial condition.
00585   
00586 }
00587 
00588   
00589 template<class Scalar> 
00590 RCP<const Thyra::ModelEvaluator<Scalar> >
00591 ForwardSensitivityStepper<Scalar>::getStateModel() const
00592 {
00593   return stateModel_;
00594 }
00595 
00596   
00597 template<class Scalar> 
00598 RCP<const ForwardSensitivityModelEvaluator<Scalar> >
00599 ForwardSensitivityStepper<Scalar>::getFwdSensModel() const
00600 {
00601   return sensModel_;
00602 }
00603 
00604   
00605 template<class Scalar> 
00606 RCP<const StateAndForwardSensitivityModelEvaluator<Scalar> >
00607 ForwardSensitivityStepper<Scalar>::getStateAndFwdSensModel() const
00608 {
00609   return stateAndSensModel_;
00610 }
00611 
00612 
00613 // Overridden from Teuchos::ParameterListAcceptor
00614 
00615 
00616 template<class Scalar> 
00617 void ForwardSensitivityStepper<Scalar>::setParameterList(
00618   RCP<Teuchos::ParameterList> const& paramList
00619   )
00620 {
00621   TEST_FOR_EXCEPT(is_null(paramList));
00622   paramList->validateParameters(*getValidParameters());
00623   this->setMyParamList(paramList);
00624   forceUpToDateW_ = paramList->get(forceUpToDateW_name_,forceUpToDateW_default_);
00625   Teuchos::readVerboseObjectSublist(&*paramList,this);
00626 }
00627 
00628 
00629 template<class Scalar> 
00630 RCP<const Teuchos::ParameterList>
00631 ForwardSensitivityStepper<Scalar>::getValidParameters() const
00632 {
00633   static RCP<const ParameterList> validPL;
00634   if (is_null(validPL) ) {
00635     RCP<ParameterList> pl = Teuchos::parameterList();
00636     pl->set( forceUpToDateW_name_, forceUpToDateW_default_,
00637       "If set to true, then the Jacobian matrix W used in the\n"
00638       "state timestep equation will be forced to be up to date\n"
00639       "with the final value for x for the nonlinear solve.  If\n"
00640       "you are willing to live with slightly less accurate sensitivities\n"
00641       "then set this to false."
00642       );
00643     Teuchos::setupVerboseObjectSublist(&*pl);
00644     validPL = pl;
00645   }
00646   return validPL;
00647 }
00648 
00649 
00650 // Overridden from StepperBase
00651 
00652 template<class Scalar> 
00653 bool ForwardSensitivityStepper<Scalar>::acceptsModel() const
00654 {
00655   return false;
00656 }
00657 
00658 template<class Scalar> 
00659 void ForwardSensitivityStepper<Scalar>::setModel(
00660   const RCP<const Thyra::ModelEvaluator<Scalar> > &model
00661   )
00662 {
00663   TEST_FOR_EXCEPT("Error, this stepper subclass does not accept a model"
00664     " as defined by the StepperBase interface!");
00665 }
00666 
00667 
00668 template<class Scalar> 
00669 RCP<const Thyra::ModelEvaluator<Scalar> >
00670 ForwardSensitivityStepper<Scalar>::getModel() const
00671 {
00672   return stateAndSensModel_;
00673 }
00674 
00675 
00676 template<class Scalar> 
00677 void ForwardSensitivityStepper<Scalar>::setInitialCondition(
00678   const Thyra::ModelEvaluatorBase::InArgs<Scalar> &state_and_sens_ic
00679   )
00680 {
00681   
00682   typedef Thyra::ModelEvaluatorBase MEB;
00683 
00684   // Get the product vectors for x_bar = [ x; s_bar ] and x_bar_dot
00685 
00686   TEST_FOR_EXCEPTION(
00687     is_null(state_and_sens_ic.get_x()), std::logic_error,
00688     "Error, the initial condition for x_bar = [ x; s_bar ] can not be null!" );
00689 
00690   const RCP<const Thyra::ProductVectorBase<Scalar> >
00691     x_bar_init = Thyra::productVectorBase<Scalar>(
00692       state_and_sens_ic.get_x()
00693       );
00694 
00695   const RCP<const Thyra::ProductVectorBase<Scalar> >
00696     x_bar_dot_init = Thyra::productVectorBase<Scalar>(
00697       state_and_sens_ic.get_x_dot()
00698       );
00699 
00700   // Remove x and x_dot from state_and_sens_ic_in to avoid cloning x and x dot!
00701   
00702   Thyra::ModelEvaluatorBase::InArgs<Scalar>
00703     state_and_sens_ic_no_x = state_and_sens_ic;
00704   state_and_sens_ic_no_x.set_x(Teuchos::null);
00705   state_and_sens_ic_no_x.set_x_dot(Teuchos::null);
00706 
00707   // Set initial condition for the state
00708 
00709   MEB::InArgs<Scalar> state_ic = stateModel_->createInArgs();
00710   state_ic.setArgs(state_and_sens_ic_no_x,true,true); // Set time, parameters etc.
00711   state_ic.set_x(x_bar_init->getVectorBlock(0)->clone_v());
00712   state_ic.set_x_dot(
00713     !is_null(x_bar_dot_init)
00714     ? x_bar_dot_init->getVectorBlock(0)->clone_v()
00715     : Teuchos::null
00716     );
00717   stateStepper_->setInitialCondition(state_ic);
00718 
00719   // Set initial condition for the sensitivities
00720   
00721   MEB::InArgs<Scalar> sens_ic = sensModel_->createInArgs();
00722   sens_ic.setArgs(state_and_sens_ic_no_x,true,true); // Set time etc.
00723   sens_ic.set_x(x_bar_init->getVectorBlock(1)->clone_v());
00724   sens_ic.set_x_dot(
00725     !is_null(x_bar_dot_init)
00726     ? x_bar_dot_init->getVectorBlock(1)->clone_v()
00727     : Teuchos::null
00728     );
00729   sensStepper_->setInitialCondition(sens_ic);
00730 
00731 }
00732  
00733 
00734 template<class Scalar> 
00735 Scalar
00736 ForwardSensitivityStepper<Scalar>::takeStep(
00737   Scalar dt, StepSizeType stepType
00738   )
00739 {
00740 
00741 #ifdef ENABLE_RYTHMOS_TIMERS
00742   TEUCHOS_FUNC_TIME_MONITOR("Rythmos:ForwardSensitivityStepper::takeStep");
00743 #endif
00744 
00745   using Teuchos::as;
00746   typedef Teuchos::ScalarTraits<Scalar> ST;
00747   typedef Teuchos::VerboseObjectTempState<InterpolationBufferBase<Scalar> > VOTSIBB;
00748   typedef Thyra::ModelEvaluatorBase MEB;
00749 
00750   RCP<Teuchos::FancyOStream> out = this->getOStream();
00751   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00752   const bool lowTrace =
00753     ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) );
00754   const bool mediumTrace =
00755     ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_MEDIUM) );
00756   Teuchos::OSTab tab(out);
00757 
00758   if (lowTrace) {
00759     *out
00760       << "\nEntering " << Teuchos::TypeNameTraits<ForwardSensitivityStepper<Scalar> >::name()
00761       << "::takeStep("<<dt<<","<<toString(stepType)<<") ...\n"; 
00762   }
00763 
00764   //
00765   // A) Compute the state timestep
00766   //
00767 
00768   if (lowTrace) {
00769     *out
00770       << "\nTaking state step using stepper : "
00771       << stateStepper_->description() << "\n";
00772   }
00773 
00774   Scalar state_dt = -1.0;
00775   {
00776 #ifdef ENABLE_RYTHMOS_TIMERS
00777     TEUCHOS_FUNC_TIME_MONITOR("Rythmos:ForwardSensitivityStepper::takeStep: stateStep");
00778 #endif
00779     VOTSIBB stateStepper_outputTempState(stateStepper_,out,verbLevel);
00780     state_dt = stateStepper_->takeStep(dt,stepType);
00781   }
00782 
00783   if (state_dt < Scalar(-ST::one())) {
00784     if (lowTrace)
00785       *out << "\nThe state stepper has failed so return a failed timestep!\n";
00786     return state_dt;
00787   }
00788   
00789   const StepStatus<Scalar> stateStepStatus = stateStepper_->getStepStatus();
00790   if (mediumTrace)
00791     *out << "\nState step status:\n" << stateStepStatus;
00792 
00793   //
00794   // B) Setup the sensitivity model for this timestep
00795   //
00796 
00797   {
00798 
00799 #ifdef ENABLE_RYTHMOS_TIMERS
00800     TEUCHOS_FUNC_TIME_MONITOR("Rythmos:ForwardSensitivityStepper::takeStep: updateSensModel");
00801 #endif
00802     
00803     TEST_FOR_EXCEPTION(
00804       stateStepStatus.stepStatus != STEP_STATUS_CONVERGED, std::logic_error,
00805       "Error, the status should be converged since a positive step size was returned!"
00806       );
00807     
00808     const Scalar curr_t = stateStepStatus.time;
00809     
00810     // Get both x and x_dot since these are needed compute other derivative
00811     // objects at these points.
00812 
00813     RCP<const Thyra::VectorBase<Scalar> > x, x_dot;
00814     get_x_and_x_dot(*stateStepper_,curr_t,&x,&x_dot);
00815     
00816     stateBasePoint_t_ = stateBasePoint_;
00817     stateBasePoint_t_.set_x_dot( x_dot );
00818     stateBasePoint_t_.set_x( x );
00819     stateBasePoint_t_.set_t( curr_t );
00820 
00821     // Grab the SingleResidualModel that was used to compute the state timestep.
00822     // From this, we can get the constants that where used to compute W!
00823     RCP<const Rythmos::SingleResidualModelEvaluatorBase<Scalar> >
00824       singleResidualModel
00825       = Teuchos::rcp_dynamic_cast<const Rythmos::SingleResidualModelEvaluatorBase<Scalar> >(
00826         stateTimeStepSolver_->getModel(), true
00827         );
00828     const Scalar
00829       coeff_x_dot = singleResidualModel->get_coeff_x_dot(),
00830       coeff_x = singleResidualModel->get_coeff_x();
00831     
00832     // Get W (and force and update if not up to date already)
00833 
00834     if (mediumTrace && forceUpToDateW_)
00835       *out << "\nForcing an update of W at the converged state timestep ...\n";
00836     
00837     RCP<Thyra::LinearOpWithSolveBase<Scalar> >
00838       W_tilde = stateTimeStepSolver_->get_nonconst_W(forceUpToDateW_);
00839     
00840     TEST_FOR_EXCEPTION(
00841       is_null(W_tilde), std::logic_error,
00842       "Error, the W from the state time step be non-null!"
00843       );
00844     
00845     sensModel_->initializeState( stateBasePoint_t_, W_tilde, coeff_x_dot, coeff_x );
00846 
00847   }
00848   // 2007/09/04: rabartl: ToDo: Move the above code that initializes the
00849   // senstivity stepper from the state stepper as a strategy object that can
00850   // be redefined.  This will allow me to make this class 100% general for all
00851   // time stepper methods ... Including explicit steppers!
00852 
00853   //
00854   // C) Compute the sensitivity timestep for the exact same timestep as was
00855   // used for the state solve.
00856   //
00857 
00858   if (lowTrace) {
00859     *out
00860       << "\nTaking sensitivity step using stepper : "
00861       << sensStepper_->description() << "\n";
00862   }
00863 
00864   Scalar sens_dt = -1.0;
00865   {
00866 
00867 #ifdef ENABLE_RYTHMOS_TIMERS
00868     TEUCHOS_FUNC_TIME_MONITOR("Rythmos:ForwardSensitivityStepper::takeStep: sensStep");
00869 #endif
00870     
00871     // Copy the step control data to make sure that the sensStepper takes the
00872     // same type of step that the statStepper took.  This is needed to ensure
00873     // that the W matrix is the same for one.
00874     sensStepper_->setStepControlData(*stateStepper_);
00875     
00876     VOTSIBB sensStepper_outputTempState(sensStepper_,out,verbLevel);
00877     
00878     sens_dt = sensStepper_->takeStep(state_dt,STEP_TYPE_FIXED);
00879     
00880   }
00881 
00882   if (mediumTrace) {
00883     const StepStatus<Scalar> sensStepStatus = sensStepper_->getStepStatus();
00884     *out
00885       << "\nSensitivity step status:\n" << sensStepStatus;
00886   }
00887   
00888   TEST_FOR_EXCEPTION(
00889     sens_dt != state_dt, std::logic_error,
00890     "Error, the sensitivity step failed for some reason.  We should\n"
00891     "just return a negative step size and reject the step but currently\n"
00892     "there is no way to roll back the state timestep it for back to\n"
00893     "the status before this function was called!"
00894     );
00895 
00896   // 2007/05/18: rabartl: ToDo: If stepType == STEP_TYPE_VARIABLE and the state
00897   // timestep sucessed but the sensitivity timestep failed, then we need to
00898   // really throw an excpetion because there is nothing that we can really do
00899   // here!
00900 
00901   // 2007/05/18: rabartl: ToDo: Replace the above std::logic_error type with
00902   // a Rythmos::CatastrophicFailure or just use Thyra::CatastrophicFailure!
00903 
00904   if (lowTrace) {
00905     *out
00906       << "\nLeaving " << Teuchos::TypeNameTraits<ForwardSensitivityStepper<Scalar> >::name()
00907       << "::takeStep("<<dt<<","<<toString(stepType)<<") ...\n"; 
00908   }
00909 
00910   return state_dt;
00911 
00912 }
00913 
00914 
00915 template<class Scalar> 
00916 const StepStatus<Scalar>
00917 ForwardSensitivityStepper<Scalar>::getStepStatus() const
00918 {
00919 
00920   const StepStatus<Scalar>
00921     stateStepStatus = stateStepper_->getStepStatus();
00922   const StepStatus<Scalar>
00923     sensStepStatus = sensStepper_->getStepStatus();
00924 
00925   StepStatus<Scalar> stepStatus;
00926   
00927   stepStatus.message = sensStepStatus.message;
00928   stepStatus.stepStatus = sensStepStatus.stepStatus;
00929   stepStatus.stepLETStatus = sensStepStatus.stepLETStatus;
00930   stepStatus.stepSize = sensStepStatus.stepSize;
00931   stepStatus.order = sensStepStatus.order;
00932   stepStatus.time = sensStepStatus.time;
00933   stepStatus.stepLETValue = sensStepStatus.stepLETValue;
00934   if (!is_null(stateStepStatus.solution) && !is_null(sensStepStatus.solution))
00935     stepStatus.solution = stateAndSensModel_->create_x_bar_vec(
00936       stateStepStatus.solution, sensStepStatus.solution
00937       );
00938   if (!is_null(stateStepStatus.solutionDot) && !is_null(sensStepStatus.solutionDot))
00939     stepStatus.solutionDot = stateAndSensModel_->create_x_bar_vec(
00940       stateStepStatus.solutionDot, sensStepStatus.solutionDot
00941       );
00942   stepStatus.extraParameters = sensStepStatus.extraParameters;
00943 
00944   return stepStatus;
00945 
00946 }
00947 
00948 
00949 // Overridden from InterpolationBufferBase
00950 
00951 
00952 template<class Scalar> 
00953 RCP<const Thyra::VectorSpaceBase<Scalar> >
00954 ForwardSensitivityStepper<Scalar>::get_x_space() const
00955 {
00956   return stateAndSensModel_->get_x_space();
00957 }
00958 
00959 
00960 template<class Scalar> 
00961 void ForwardSensitivityStepper<Scalar>::addPoints(
00962   const Array<Scalar>& time_vec,
00963   const Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >& x_vec,
00964   const Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >& xdot_vec
00965   )
00966 {
00967   TEST_FOR_EXCEPT("Not implemented addPoints(...) yet but we could if we wanted!");
00968 }
00969 
00970 
00971 template<class Scalar> 
00972 TimeRange<Scalar>
00973 ForwardSensitivityStepper<Scalar>::getTimeRange() const
00974 {
00975   return sensStepper_->getTimeRange();
00976 }
00977 
00978 
00979 template<class Scalar> 
00980 void ForwardSensitivityStepper<Scalar>::getPoints(
00981   const Array<Scalar>& time_vec,
00982   Array<RCP<const Thyra::VectorBase<Scalar> > >* x_bar_vec,
00983   Array<RCP<const Thyra::VectorBase<Scalar> > >* x_bar_dot_vec,
00984   Array<ScalarMag>* accuracy_vec
00985   ) const
00986 {
00987 
00988   using Teuchos::as;
00989 
00990 #ifdef TEUCHOS_DEBUG
00991   TEST_FOR_EXCEPT( as<int>(time_vec.size()) == 0 );
00992 #endif
00993 
00994   const int numTimePoints = time_vec.size();
00995 
00996   if (x_bar_vec)
00997     x_bar_vec->clear();
00998 
00999   if (x_bar_dot_vec)
01000     x_bar_dot_vec->clear();
01001   
01002   Array<RCP<const Thyra::VectorBase<Scalar> > >
01003     x_vec, x_dot_vec;
01004 
01005   stateStepper_->getPoints(
01006     time_vec,
01007     x_bar_vec ? &x_vec: 0,
01008     x_bar_dot_vec ? &x_dot_vec: 0,
01009     0 // Ignoring accuracy from state for now!
01010     );
01011   
01012   Array<RCP<const Thyra::VectorBase<Scalar> > >
01013     s_bar_vec, s_bar_dot_vec;
01014 
01015   sensStepper_->getPoints(
01016     time_vec,
01017     x_bar_vec ? &s_bar_vec: 0,
01018     x_bar_dot_vec ? &s_bar_dot_vec: 0,
01019     accuracy_vec
01020     );
01021   
01022   if ( x_bar_vec ) {
01023     for ( int i = 0; i < numTimePoints; ++i ) {
01024       x_bar_vec->push_back(
01025         stateAndSensModel_->create_x_bar_vec(x_vec[i],s_bar_vec[i])
01026         );
01027     }
01028   }
01029   
01030   if ( x_bar_dot_vec ) {
01031     for ( int i = 0; i < numTimePoints; ++i ) {
01032       x_bar_dot_vec->push_back(
01033         stateAndSensModel_->create_x_bar_vec(x_dot_vec[i],s_bar_dot_vec[i])
01034         );
01035     }
01036   }
01037 
01038 }
01039 
01040 
01041 template<class Scalar>
01042 void ForwardSensitivityStepper<Scalar>::getNodes(
01043   Array<Scalar>* time_vec
01044   ) const
01045 {
01046   TEST_FOR_EXCEPT("Not implemented yet but we can!");
01047 }
01048 
01049 
01050 template<class Scalar> 
01051 void ForwardSensitivityStepper<Scalar>::removeNodes(
01052   Array<Scalar>& time_vec
01053   )
01054 {
01055   TEST_FOR_EXCEPT("Not implemented yet but we can!");
01056 }
01057 
01058 
01059 template<class Scalar> 
01060 int ForwardSensitivityStepper<Scalar>::getOrder() const
01061 {
01062   return sensStepper_->getOrder();
01063   // Note: This assumes that stateStepper will have the same order!
01064 }
01065 
01066 
01067 } // namespace Rythmos
01068 
01069 
01070 #endif //RYTHMOS_FORWARD_SENSITIVITY_STEPPER_HPP

Generated on Tue Oct 20 12:46:08 2009 for Rythmos - Transient Integration for Differential Equations by doxygen 1.4.7