Rythmos - Transient Integration for Differential Equations Version of the Day
Rythmos_DiagonalImplicitRKModelEvaluator.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., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
00023 // USA
00024 // Questions? Contact Todd S. Coffey (tscoffe@sandia.gov)
00025 //
00026 // ***********************************************************************
00027 //@HEADER
00028 
00029 
00030 #ifndef RYTHMOS_DIAGONALIMPLICITRK_MODEL_EVALUATOR_HPP
00031 #define RYTHMOS_DIAGONALIMPLICITRK_MODEL_EVALUATOR_HPP
00032 
00033 
00034 #include "Rythmos_Types.hpp"
00035 #include "Rythmos_RKButcherTableauHelpers.hpp"
00036 #include "Thyra_StateFuncModelEvaluatorBase.hpp"
00037 #include "Thyra_ModelEvaluatorHelpers.hpp"
00038 #include "Thyra_ModelEvaluatorDelegatorBase.hpp"
00039 #include "Thyra_DefaultProductVectorSpace.hpp"
00040 #include "Thyra_DefaultBlockedLinearOp.hpp"
00041 #include "Thyra_VectorStdOps.hpp"
00042 #include "Thyra_TestingTools.hpp"
00043 #include "Teuchos_as.hpp"
00044 #include "Teuchos_SerialDenseMatrix.hpp"
00045 #include "Teuchos_SerialDenseVector.hpp"
00046 #include "Teuchos_Assert.hpp"
00047 
00048 
00049 namespace Rythmos {
00050 
00051 
00053 template<class Scalar>
00054 class DiagonalImplicitRKModelEvaluator : virtual public Thyra::StateFuncModelEvaluatorBase<Scalar>
00055 {
00056 public:
00057 
00060 
00062   DiagonalImplicitRKModelEvaluator();
00063 
00065   void initializeDIRKModel(
00066     const RCP<const Thyra::ModelEvaluator<Scalar> >& daeModel,
00067     const Thyra::ModelEvaluatorBase::InArgs<Scalar>& basePoint,
00068     const RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> >& dirk_W_factory,
00069     const RCP<const RKButcherTableauBase<Scalar> >& irkButcherTableau
00070     );
00071 
00073   void setTimeStepPoint(
00074     const RCP<const Thyra::VectorBase<Scalar> > &x_old,
00075     const Scalar &t_old,
00076     const Scalar &delta_t
00077     );
00078 
00080   void setStageSolution(
00081       int stageNumber,
00082       const Thyra::VectorBase<Scalar>& stage_solution
00083       );
00084 
00086   void setCurrentStage(
00087       int currentStage
00088       );
00089 
00091 
00094 
00096   RCP<const Thyra::VectorSpaceBase<Scalar> > get_x_space() const;
00098   RCP<const Thyra::VectorSpaceBase<Scalar> > get_f_space() const;
00100   RCP<Thyra::LinearOpBase<Scalar> > create_W_op() const;
00102   RCP<const Thyra::LinearOpWithSolveFactoryBase<Scalar> > get_W_factory() const;
00104   Thyra::ModelEvaluatorBase::InArgs<Scalar> getNominalValues() const;
00106   Thyra::ModelEvaluatorBase::InArgs<Scalar> createInArgs() const;
00107 
00109 
00110 private:
00111 
00114 
00116   Thyra::ModelEvaluatorBase::OutArgs<Scalar> createOutArgsImpl() const;
00118   void evalModelImpl(
00119     const Thyra::ModelEvaluatorBase::InArgs<Scalar>& inArgs,
00120     const Thyra::ModelEvaluatorBase::OutArgs<Scalar>& outArgs
00121     ) const;
00122 
00124 
00125 
00126 private:
00127 
00128   RCP<const Thyra::ModelEvaluator<Scalar> > daeModel_;
00129   Thyra::ModelEvaluatorBase::InArgs<Scalar> basePoint_;
00130   RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> > dirk_W_factory_;
00131   RCP<const RKButcherTableauBase<Scalar> > dirkButcherTableau_;
00132 
00133   Thyra::ModelEvaluatorBase::InArgs<Scalar> inArgs_;
00134   Thyra::ModelEvaluatorBase::OutArgs<Scalar> outArgs_;
00135   Thyra::ModelEvaluatorBase::InArgs<Scalar> nominalValues_;
00136 
00137   RCP<const Thyra::ProductVectorSpaceBase<Scalar> > stage_space_;
00138   RCP<Thyra::ProductVectorBase<Scalar> > stage_derivatives_;
00139 
00140   bool setTimeStepPointCalled_;
00141   RCP<const Thyra::VectorBase<Scalar> > x_old_;
00142   Scalar t_old_;
00143   Scalar delta_t_;
00144   int currentStage_;
00145 
00146   bool isInitialized_;
00147 
00148 };
00149 
00150 
00155 template<class Scalar>
00156 RCP<DiagonalImplicitRKModelEvaluator<Scalar> >
00157 diagonalImplicitRKModelEvaluator(
00158   const RCP<const Thyra::ModelEvaluator<Scalar> >& daeModel,
00159   const Thyra::ModelEvaluatorBase::InArgs<Scalar> &basePoint,
00160   const RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> > &dirk_W_factory,
00161   const RCP<const RKButcherTableauBase<Scalar> > &irkButcherTableau
00162   )
00163 {
00164   RCP<DiagonalImplicitRKModelEvaluator<Scalar> >
00165     dirkModel = rcp(new DiagonalImplicitRKModelEvaluator<Scalar>());
00166   dirkModel->initializeDIRKModel(daeModel,basePoint,dirk_W_factory,irkButcherTableau);
00167   return dirkModel;
00168 }
00169 
00170 
00171 // ///////////////////////
00172 // Definition
00173 
00174 
00175 // Constructors/initializers/accessors
00176 
00177 
00178 template<class Scalar>
00179 DiagonalImplicitRKModelEvaluator<Scalar>::DiagonalImplicitRKModelEvaluator()
00180 {
00181   setTimeStepPointCalled_ = false;
00182   isInitialized_ = false;
00183   currentStage_ = -1;
00184 }
00185 
00186 
00187 // Overridden from ImplicitRKModelEvaluatorBase
00188 
00189 
00190 template<class Scalar>
00191 void DiagonalImplicitRKModelEvaluator<Scalar>::initializeDIRKModel(
00192   const RCP<const Thyra::ModelEvaluator<Scalar> >& daeModel,
00193   const Thyra::ModelEvaluatorBase::InArgs<Scalar>& basePoint,
00194   const RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> >& dirk_W_factory,
00195   const RCP<const RKButcherTableauBase<Scalar> >& irkButcherTableau
00196   )
00197 {
00198   typedef ScalarTraits<Scalar> ST;
00199   // ToDo: Assert input arguments!
00200   // How do I verify the basePoint is an authentic InArgs from daeModel?
00201   TEUCHOS_TEST_FOR_EXCEPTION( 
00202       is_null(basePoint.get_x()), 
00203       std::logic_error,
00204       "Error!  The basepoint x vector is null!"
00205       );
00206   TEUCHOS_TEST_FOR_EXCEPTION( 
00207       is_null(daeModel), 
00208       std::logic_error,
00209       "Error!  The model evaluator pointer is null!"
00210       );
00211   TEUCHOS_TEST_FOR_EXCEPTION( 
00212       !daeModel->get_x_space()->isCompatible(*(basePoint.get_x()->space())), 
00213       std::logic_error,
00214       "Error!  The basepoint input arguments are incompatible with the model evaluator vector space!"
00215       );
00216   //TEUCHOS_TEST_FOR_EXCEPT(is_null(dirk_W_factory));
00217 
00218   daeModel_ = daeModel;
00219   basePoint_ = basePoint;
00220   dirk_W_factory_ = dirk_W_factory;
00221   dirkButcherTableau_ = irkButcherTableau;
00222 
00223   const int numStages = dirkButcherTableau_->numStages();
00224 
00225   using Teuchos::rcp_dynamic_cast;
00226   stage_space_ = productVectorSpace(daeModel_->get_x_space(),numStages);
00227   RCP<const Thyra::VectorSpaceBase<Scalar> > vs = rcp_dynamic_cast<const Thyra::VectorSpaceBase<Scalar> >(stage_space_,true);
00228   stage_derivatives_ = rcp_dynamic_cast<Thyra::ProductVectorBase<Scalar> >(createMember(vs),true);
00229   Thyra::V_S( rcp_dynamic_cast<Thyra::VectorBase<Scalar> >(stage_derivatives_).ptr(),ST::zero());
00230 
00231   // Set up prototypical InArgs
00232   {
00233     typedef Thyra::ModelEvaluatorBase MEB;
00234     MEB::InArgsSetup<Scalar> inArgs;
00235     inArgs.setModelEvalDescription(this->description());
00236     inArgs.setSupports(MEB::IN_ARG_x);
00237     inArgs_ = inArgs;
00238   }
00239   // Set up prototypical OutArgs
00240   {
00241     typedef Thyra::ModelEvaluatorBase MEB;
00242     MEB::OutArgsSetup<Scalar> outArgs;
00243     outArgs.setModelEvalDescription(this->description());
00244     outArgs.setSupports(MEB::OUT_ARG_f);
00245     outArgs.setSupports(MEB::OUT_ARG_W_op);
00246     outArgs_ = outArgs;
00247   }
00248   // Set up nominal values
00249   nominalValues_ = inArgs_;
00250 
00251   isInitialized_ = true;
00252 }
00253 
00254 
00255 template<class Scalar>
00256 void DiagonalImplicitRKModelEvaluator<Scalar>::setTimeStepPoint(
00257   const RCP<const Thyra::VectorBase<Scalar> > &x_old,
00258   const Scalar &t_old,
00259   const Scalar &delta_t
00260   )
00261 {
00262   typedef ScalarTraits<Scalar> ST;
00263   TEUCHOS_TEST_FOR_EXCEPT( is_null(x_old) );
00264   TEUCHOS_TEST_FOR_EXCEPT( delta_t <= ST::zero() );
00265   // Verify x_old is compatible with the vector space in the DAE Model.
00266   TEUCHOS_TEST_FOR_EXCEPTION(!daeModel_->get_x_space()->isCompatible(*(x_old->space())), std::logic_error,
00267       "Error!  The incoming VectorBase object is not compatible with the DAE model that was provided at initialization!");
00268   x_old_ = x_old;
00269   t_old_ = t_old;
00270   delta_t_ = delta_t;
00271   setTimeStepPointCalled_ = true;
00272 }
00273 
00274 template<class Scalar>
00275 void DiagonalImplicitRKModelEvaluator<Scalar>::setStageSolution(
00276       int stageNumber,
00277       const Thyra::VectorBase<Scalar>& stage_solution
00278       )
00279 {
00280   TEUCHOS_TEST_FOR_EXCEPT(stageNumber < 0);
00281   TEUCHOS_TEST_FOR_EXCEPT(stageNumber >= dirkButcherTableau_->numStages());
00282   Thyra::V_V(stage_derivatives_->getNonconstVectorBlock(stageNumber).ptr(),stage_solution);
00283 }
00284 
00285 template<class Scalar>
00286 void DiagonalImplicitRKModelEvaluator<Scalar>::setCurrentStage(
00287       int currentStage
00288       )
00289 {
00290   TEUCHOS_TEST_FOR_EXCEPT(currentStage < 0);
00291   TEUCHOS_TEST_FOR_EXCEPT(currentStage >= dirkButcherTableau_->numStages());
00292   currentStage_ = currentStage;
00293 }
00294 
00295 
00296 
00297 // Overridden from ModelEvaluator
00298 
00299 
00300 template<class Scalar>
00301 RCP<const Thyra::VectorSpaceBase<Scalar> >
00302 DiagonalImplicitRKModelEvaluator<Scalar>::get_x_space() const
00303 {
00304   TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
00305       "Error, initializeDIRKModel must be called first!\n"
00306       );
00307   return daeModel_->get_x_space();
00308 }
00309 
00310 
00311 template<class Scalar>
00312 RCP<const Thyra::VectorSpaceBase<Scalar> >
00313 DiagonalImplicitRKModelEvaluator<Scalar>::get_f_space() const
00314 {
00315   TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
00316       "Error, initializeDIRKModel must be called first!\n"
00317       );
00318   return daeModel_->get_f_space();
00319 }
00320 
00321 
00322 template<class Scalar>
00323 RCP<Thyra::LinearOpBase<Scalar> >
00324 DiagonalImplicitRKModelEvaluator<Scalar>::create_W_op() const
00325 {
00326   TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
00327       "Error, initializeDIRKModel must be called first!\n"
00328       );
00329   return daeModel_->create_W_op();
00330 }
00331 
00332 
00333 template<class Scalar>
00334 RCP<const Thyra::LinearOpWithSolveFactoryBase<Scalar> >
00335 DiagonalImplicitRKModelEvaluator<Scalar>::get_W_factory() const
00336 {
00337   TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
00338       "Error, initializeDIRKModel must be called first!\n"
00339       );
00340   return daeModel_->get_W_factory();
00341 }  
00342 
00343 
00344 template<class Scalar>
00345 Thyra::ModelEvaluatorBase::InArgs<Scalar>
00346 DiagonalImplicitRKModelEvaluator<Scalar>::getNominalValues() const
00347 {
00348   TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
00349       "Error, initializeDIRKModel must be called first!\n"
00350       );
00351   return nominalValues_;
00352 }
00353 
00354 
00355 template<class Scalar>
00356 Thyra::ModelEvaluatorBase::InArgs<Scalar>
00357 DiagonalImplicitRKModelEvaluator<Scalar>::createInArgs() const
00358 {
00359   TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
00360       "Error, initializeDIRKModel must be called first!\n"
00361       );
00362   return inArgs_;
00363 }
00364 
00365 
00366 // Private functions overridden from ModelEvaluatorDefaultBase
00367 
00368 
00369 template<class Scalar>
00370 Thyra::ModelEvaluatorBase::OutArgs<Scalar>
00371 DiagonalImplicitRKModelEvaluator<Scalar>::createOutArgsImpl() const
00372 {
00373   TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
00374       "Error, initializeDIRKModel must be called first!\n"
00375       );
00376   return outArgs_;
00377 }
00378 
00379 
00380 template<class Scalar>
00381 void DiagonalImplicitRKModelEvaluator<Scalar>::evalModelImpl(
00382   const Thyra::ModelEvaluatorBase::InArgs<Scalar>& inArgs_stage,
00383   const Thyra::ModelEvaluatorBase::OutArgs<Scalar>& outArgs_stage
00384   ) const
00385 {
00386 
00387   typedef ScalarTraits<Scalar> ST;
00388   typedef Thyra::ModelEvaluatorBase MEB;
00389 
00390   TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
00391       "Error!  initializeDIRKModel must be called before evalModel\n"
00392       );
00393 
00394   TEUCHOS_TEST_FOR_EXCEPTION( !setTimeStepPointCalled_, std::logic_error,
00395       "Error!  setTimeStepPoint must be called before evalModel"
00396       );
00397 
00398   TEUCHOS_TEST_FOR_EXCEPTION( currentStage_ == -1, std::logic_error,
00399       "Error!  setCurrentStage must be called before evalModel"
00400       );
00401 
00402   THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_GEN_BEGIN(
00403     "Rythmos::DiagonalImplicitRKModelEvaluator",inArgs_stage,outArgs_stage,daeModel_
00404     );
00405 
00406   //
00407   // A) Unwrap the inArgs and outArgs 
00408   //
00409 
00410   const RCP<const Thyra::VectorBase<Scalar> > x_in = inArgs_stage.get_x();
00411   const RCP<Thyra::VectorBase<Scalar> > f_out = outArgs_stage.get_f();
00412   const RCP<Thyra::LinearOpBase<Scalar> > W_op_out = outArgs_stage.get_W_op();
00413 
00414   //
00415   // B) Assemble f_out and W_op_out for given stage
00416   //
00417 
00418   MEB::InArgs<Scalar> daeInArgs = daeModel_->createInArgs();
00419   MEB::OutArgs<Scalar> daeOutArgs = daeModel_->createOutArgs();
00420   const RCP<Thyra::VectorBase<Scalar> > x_i = createMember(daeModel_->get_x_space());
00421   daeInArgs.setArgs(basePoint_);
00422   
00423   // B.1) Setup the DAE's inArgs for stage f(currentStage_) ...
00424   V_V(stage_derivatives_->getNonconstVectorBlock(currentStage_).ptr(),*x_in);
00425   assembleIRKState( currentStage_, dirkButcherTableau_->A(), delta_t_, *x_old_, *stage_derivatives_, outArg(*x_i) );
00426   daeInArgs.set_x( x_i );
00427   daeInArgs.set_x_dot( x_in );
00428   daeInArgs.set_t( t_old_ + dirkButcherTableau_->c()(currentStage_) * delta_t_ );
00429   daeInArgs.set_alpha(ST::one());
00430   daeInArgs.set_beta( delta_t_ * dirkButcherTableau_->A()(currentStage_,currentStage_) );
00431 
00432   // B.2) Setup the DAE's outArgs for stage f(i) ...
00433   if (!is_null(f_out))
00434     daeOutArgs.set_f( f_out );
00435   if (!is_null(W_op_out))
00436     daeOutArgs.set_W_op(W_op_out);
00437 
00438   // B.3) Compute f_out(i) and/or W_op_out ...
00439   daeModel_->evalModel( daeInArgs, daeOutArgs );
00440   daeOutArgs.set_f(Teuchos::null);
00441   daeOutArgs.set_W_op(Teuchos::null);
00442   
00443   THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_END();
00444   
00445 }
00446 
00447 
00448 } // namespace Rythmos
00449 
00450 
00451 #endif // RYTHMOS_DIAGONALIMPLICITRK_MODEL_EVALUATOR_HPP
 All Classes Functions Variables Typedefs Friends