Rythmos - Transient Integration for Differential Equations Version of the Day
Rythmos_ImplicitRKModelEvaluator.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_IMPLICITRK_MODEL_EVALUATOR_HPP
00031 #define RYTHMOS_IMPLICITRK_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 ImplicitRKModelEvaluator : virtual public Thyra::StateFuncModelEvaluatorBase<Scalar>
00055 {
00056 public:
00057 
00060 
00062   ImplicitRKModelEvaluator();
00063 
00065   void initializeIRKModel(
00066     const RCP<const Thyra::ModelEvaluator<Scalar> >& daeModel,
00067     const Thyra::ModelEvaluatorBase::InArgs<Scalar>& basePoint,
00068     const RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> >& irk_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 
00083 
00085   RCP<const Thyra::VectorSpaceBase<Scalar> > get_x_space() const;
00087   RCP<const Thyra::VectorSpaceBase<Scalar> > get_f_space() const;
00089   RCP<Thyra::LinearOpBase<Scalar> > create_W_op() const;
00091   RCP<const Thyra::LinearOpWithSolveFactoryBase<Scalar> > get_W_factory() const;
00093   Thyra::ModelEvaluatorBase::InArgs<Scalar> getNominalValues() const;
00095   Thyra::ModelEvaluatorBase::InArgs<Scalar> createInArgs() const;
00096 
00098 
00099 private:
00100 
00103 
00105   Thyra::ModelEvaluatorBase::OutArgs<Scalar> createOutArgsImpl() const;
00107   void evalModelImpl(
00108     const Thyra::ModelEvaluatorBase::InArgs<Scalar>& inArgs,
00109     const Thyra::ModelEvaluatorBase::OutArgs<Scalar>& outArgs
00110     ) const;
00111 
00113 
00114 
00115 private:
00116 
00117   RCP<const Thyra::ModelEvaluator<Scalar> > daeModel_;
00118   Thyra::ModelEvaluatorBase::InArgs<Scalar> basePoint_;
00119   RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> > irk_W_factory_;
00120   RCP<const RKButcherTableauBase<Scalar> > irkButcherTableau_;
00121 
00122   Thyra::ModelEvaluatorBase::InArgs<Scalar> inArgs_;
00123   Thyra::ModelEvaluatorBase::OutArgs<Scalar> outArgs_;
00124   Thyra::ModelEvaluatorBase::InArgs<Scalar> nominalValues_;
00125 
00126   RCP<const Thyra::ProductVectorSpaceBase<Scalar> > x_bar_space_;
00127   RCP<const Thyra::ProductVectorSpaceBase<Scalar> > f_bar_space_;
00128   RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> > W_bar_factory_;
00129 
00130   bool setTimeStepPointCalled_;
00131   RCP<const Thyra::VectorBase<Scalar> > x_old_;
00132   Scalar t_old_;
00133   Scalar delta_t_;
00134 
00135   bool isInitialized_;
00136 
00137 };
00138 
00139 
00144 template<class Scalar>
00145 RCP<ImplicitRKModelEvaluator<Scalar> >
00146 implicitRKModelEvaluator(
00147   const RCP<const Thyra::ModelEvaluator<Scalar> >& daeModel,
00148   const Thyra::ModelEvaluatorBase::InArgs<Scalar>& basePoint,
00149   const RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> >& irk_W_factory,
00150   const RCP<const RKButcherTableauBase<Scalar> >& irkButcherTableau
00151   )
00152 {
00153   RCP<ImplicitRKModelEvaluator<Scalar> >
00154     irkModel = rcp(new ImplicitRKModelEvaluator<Scalar>());
00155   irkModel->initializeIRKModel(daeModel,basePoint,irk_W_factory,irkButcherTableau);
00156   return irkModel;
00157 }
00158 
00159 
00160 // ///////////////////////
00161 // Definition
00162 
00163 
00164 // Constructors/initializers/accessors
00165 
00166 
00167 template<class Scalar>
00168 ImplicitRKModelEvaluator<Scalar>::ImplicitRKModelEvaluator()
00169 {
00170   setTimeStepPointCalled_ = false;
00171   isInitialized_ = false;
00172 }
00173 
00174 
00175 // Overridden from ImplicitRKModelEvaluatorBase
00176 
00177 
00178 template<class Scalar>
00179 void ImplicitRKModelEvaluator<Scalar>::initializeIRKModel(
00180   const RCP<const Thyra::ModelEvaluator<Scalar> >& daeModel,
00181   const Thyra::ModelEvaluatorBase::InArgs<Scalar>& basePoint,
00182   const RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> >& irk_W_factory,
00183   const RCP<const RKButcherTableauBase<Scalar> >& irkButcherTableau
00184   )
00185 {
00186   // ToDo: Assert input arguments!
00187   // How do I verify the basePoint is an authentic InArgs from daeModel?
00188   TEUCHOS_TEST_FOR_EXCEPTION( 
00189       is_null(basePoint.get_x()), 
00190       std::logic_error,
00191       "Error!  The basepoint x vector is null!"
00192       );
00193   TEUCHOS_TEST_FOR_EXCEPTION( 
00194       is_null(daeModel), 
00195       std::logic_error,
00196       "Error!  The model evaluator pointer is null!"
00197       );
00198   TEUCHOS_TEST_FOR_EXCEPTION( 
00199       !daeModel->get_x_space()->isCompatible(*(basePoint.get_x()->space())), 
00200       std::logic_error,
00201       "Error!  The basepoint input arguments are incompatible with the model evaluator vector space!"
00202       );
00203   TEUCHOS_TEST_FOR_EXCEPT(is_null(irk_W_factory));
00204 
00205   daeModel_ = daeModel;
00206   basePoint_ = basePoint;
00207   irk_W_factory_ = irk_W_factory;
00208   irkButcherTableau_ = irkButcherTableau;
00209 
00210   const int numStages = irkButcherTableau_->numStages();
00211 
00212   x_bar_space_ = productVectorSpace(daeModel_->get_x_space(),numStages);
00213   f_bar_space_ = productVectorSpace(daeModel_->get_f_space(),numStages);
00214 
00215   // HACK! Remove the preconditioner factory for now!
00216   if (irk_W_factory_->acceptsPreconditionerFactory())
00217     irk_W_factory_->unsetPreconditionerFactory();
00218 
00219   // ToDo: create the block diagonal preconditioner factory and set this on
00220   // irk_W_factory_!
00221   
00222   // Set up prototypical InArgs
00223   {
00224     typedef Thyra::ModelEvaluatorBase MEB;
00225     MEB::InArgsSetup<Scalar> inArgs;
00226     inArgs.setModelEvalDescription(this->description());
00227     inArgs.setSupports(MEB::IN_ARG_x);
00228     inArgs_ = inArgs;
00229   }
00230   // Set up prototypical OutArgs
00231   {
00232     typedef Thyra::ModelEvaluatorBase MEB;
00233     MEB::OutArgsSetup<Scalar> outArgs;
00234     outArgs.setModelEvalDescription(this->description());
00235     outArgs.setSupports(MEB::OUT_ARG_f);
00236     outArgs.setSupports(MEB::OUT_ARG_W_op);
00237     outArgs_ = outArgs;
00238   }
00239   // Set up nominal values
00240   nominalValues_ = inArgs_;
00241 
00242   isInitialized_ = true;
00243 }
00244 
00245 
00246 template<class Scalar>
00247 void ImplicitRKModelEvaluator<Scalar>::setTimeStepPoint(
00248   const RCP<const Thyra::VectorBase<Scalar> > &x_old,
00249   const Scalar &t_old,
00250   const Scalar &delta_t
00251   )
00252 {
00253   typedef ScalarTraits<Scalar> ST;
00254   TEUCHOS_TEST_FOR_EXCEPT( is_null(x_old) );
00255   TEUCHOS_TEST_FOR_EXCEPT( delta_t <= ST::zero() );
00256   // Verify x_old is compatible with the vector space in the DAE Model.
00257   TEUCHOS_TEST_FOR_EXCEPTION(!daeModel_->get_x_space()->isCompatible(*(x_old->space())), std::logic_error,
00258       "Error!  The incoming VectorBase object is not compatible with the DAE model that was provided at initialization!");
00259   x_old_ = x_old;
00260   t_old_ = t_old;
00261   delta_t_ = delta_t;
00262   setTimeStepPointCalled_ = true;
00263 }
00264 
00265 
00266 // Overridden from ModelEvaluator
00267 
00268 
00269 template<class Scalar>
00270 RCP<const Thyra::VectorSpaceBase<Scalar> >
00271 ImplicitRKModelEvaluator<Scalar>::get_x_space() const
00272 {
00273   TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
00274       "Error, initializeIRKModel must be called first!\n"
00275       );
00276   return x_bar_space_;
00277 }
00278 
00279 
00280 template<class Scalar>
00281 RCP<const Thyra::VectorSpaceBase<Scalar> >
00282 ImplicitRKModelEvaluator<Scalar>::get_f_space() const
00283 {
00284   TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
00285       "Error, initializeIRKModel must be called first!\n"
00286       );
00287   return f_bar_space_;
00288 }
00289 
00290 
00291 template<class Scalar>
00292 RCP<Thyra::LinearOpBase<Scalar> >
00293 ImplicitRKModelEvaluator<Scalar>::create_W_op() const
00294 {
00295   TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
00296       "Error, initializeIRKModel must be called first!\n"
00297       );
00298   // Create the block structure for W_op_bar right away!
00299   const int numStages = irkButcherTableau_->numStages();
00300   RCP<Thyra::PhysicallyBlockedLinearOpBase<Scalar> >
00301     W_op_bar = Thyra::defaultBlockedLinearOp<Scalar>();
00302   W_op_bar->beginBlockFill( f_bar_space_, x_bar_space_ );
00303   for ( int i = 0; i < numStages; ++i )
00304     for ( int j = 0; j < numStages; ++j )
00305       W_op_bar->setNonconstBlock( i, j, daeModel_->create_W_op() );
00306   W_op_bar->endBlockFill();
00307   return W_op_bar;
00308 }
00309 
00310 
00311 template<class Scalar>
00312 RCP<const Thyra::LinearOpWithSolveFactoryBase<Scalar> >
00313 ImplicitRKModelEvaluator<Scalar>::get_W_factory() const
00314 {
00315   TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
00316       "Error, initializeIRKModel must be called first!\n"
00317       );
00318   return irk_W_factory_;
00319 }  
00320 
00321 
00322 template<class Scalar>
00323 Thyra::ModelEvaluatorBase::InArgs<Scalar>
00324 ImplicitRKModelEvaluator<Scalar>::getNominalValues() const
00325 {
00326   TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
00327       "Error, initializeIRKModel must be called first!\n"
00328       );
00329   return nominalValues_;
00330 }
00331 
00332 
00333 template<class Scalar>
00334 Thyra::ModelEvaluatorBase::InArgs<Scalar>
00335 ImplicitRKModelEvaluator<Scalar>::createInArgs() const
00336 {
00337   TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
00338       "Error, initializeIRKModel must be called first!\n"
00339       );
00340   return inArgs_;
00341 }
00342 
00343 
00344 // Private functions overridden from ModelEvaluatorDefaultBase
00345 
00346 
00347 template<class Scalar>
00348 Thyra::ModelEvaluatorBase::OutArgs<Scalar>
00349 ImplicitRKModelEvaluator<Scalar>::createOutArgsImpl() const
00350 {
00351   TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
00352       "Error, initializeIRKModel must be called first!\n"
00353       );
00354   return outArgs_;
00355 }
00356 
00357 
00358 template<class Scalar>
00359 void ImplicitRKModelEvaluator<Scalar>::evalModelImpl(
00360   const Thyra::ModelEvaluatorBase::InArgs<Scalar>& inArgs_bar,
00361   const Thyra::ModelEvaluatorBase::OutArgs<Scalar>& outArgs_bar
00362   ) const
00363 {
00364 
00365   using Teuchos::rcp_dynamic_cast;
00366   typedef ScalarTraits<Scalar> ST;
00367   typedef Thyra::ModelEvaluatorBase MEB;
00368   typedef Thyra::VectorBase<Scalar> VB;
00369   typedef Thyra::ProductVectorBase<Scalar> PVB;
00370   typedef Thyra::BlockedLinearOpBase<Scalar> BLWB;
00371 
00372   TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
00373       "Error!  initializeIRKModel must be called before evalModel\n"
00374       );
00375 
00376   TEUCHOS_TEST_FOR_EXCEPTION( !setTimeStepPointCalled_, std::logic_error,
00377       "Error!  setTimeStepPoint must be called before evalModel"
00378       );
00379 
00380   THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_GEN_BEGIN(
00381     "Rythmos::ImplicitRKModelEvaluator",inArgs_bar,outArgs_bar,daeModel_
00382     );
00383 
00384   //
00385   // A) Unwrap the inArgs and outArgs to get at product vectors and block op
00386   //
00387 
00388   const RCP<const PVB> x_bar = rcp_dynamic_cast<const PVB>(inArgs_bar.get_x(), true);
00389   const RCP<PVB> f_bar = rcp_dynamic_cast<PVB>(outArgs_bar.get_f(), true);
00390   const RCP<BLWB> W_op_bar = rcp_dynamic_cast<BLWB>(outArgs_bar.get_W_op(), true);
00391 
00392   //
00393   // B) Assemble f_bar and W_op_bar by looping over stages
00394   //
00395 
00396   MEB::InArgs<Scalar> daeInArgs = daeModel_->createInArgs();
00397   MEB::OutArgs<Scalar> daeOutArgs = daeModel_->createOutArgs();
00398   const RCP<VB> x_i = createMember(daeModel_->get_x_space());
00399   daeInArgs.setArgs(basePoint_);
00400   
00401   const int numStages = irkButcherTableau_->numStages();
00402 
00403   for ( int i = 0; i < numStages; ++i ) {
00404 
00405     // B.1) Setup the DAE's inArgs for stage f(i) ...
00406     assembleIRKState( i, irkButcherTableau_->A(), delta_t_, *x_old_, *x_bar, outArg(*x_i) );
00407     daeInArgs.set_x( x_i );
00408     daeInArgs.set_x_dot( x_bar->getVectorBlock(i) );
00409     daeInArgs.set_t( t_old_ + irkButcherTableau_->c()(i) * delta_t_ );
00410     Scalar alpha = ST::zero();
00411     if (i == 0) {
00412       alpha = ST::one();
00413     } else {
00414       alpha = ST::zero();
00415     }
00416     Scalar beta = delta_t_ * irkButcherTableau_->A()(i,0);
00417     daeInArgs.set_alpha( alpha );
00418     daeInArgs.set_beta( beta );
00419 
00420     // B.2) Setup the DAE's outArgs for stage f(i) ...
00421     if (!is_null(f_bar))
00422       daeOutArgs.set_f( f_bar->getNonconstVectorBlock(i) );
00423     if (!is_null(W_op_bar)) {
00424       daeOutArgs.set_W_op(W_op_bar->getNonconstBlock(i,0));
00425     }
00426 
00427     // B.3) Compute f_bar(i) and/or W_op_bar(i,0) ...
00428     daeModel_->evalModel( daeInArgs, daeOutArgs );
00429     daeOutArgs.set_f(Teuchos::null);
00430     daeOutArgs.set_W_op(Teuchos::null);
00431     
00432     // B.4) Evaluate the rest of the W_op_bar(i,j=1...numStages-1) ...
00433     if (!is_null(W_op_bar)) {
00434       for ( int j = 1; j < numStages; ++j ) {
00435         alpha = ST::zero();
00436         if (i == j) {
00437           alpha = ST::one();
00438         } else {
00439           alpha = ST::zero();
00440         }
00441         beta = delta_t_ * irkButcherTableau_->A()(i,j);
00442         daeInArgs.set_alpha( alpha );
00443         daeInArgs.set_beta( beta );
00444         daeOutArgs.set_W_op(W_op_bar->getNonconstBlock(i,j));
00445         daeModel_->evalModel( daeInArgs, daeOutArgs );
00446         daeOutArgs.set_W_op(Teuchos::null);
00447       }
00448     }
00449 
00450   }
00451   
00452   THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_END();
00453   
00454 }
00455 
00456 
00457 } // namespace Rythmos
00458 
00459 
00460 #endif // RYTHMOS_IMPLICITRK_MODEL_EVALUATOR_HPP
 All Classes Functions Variables Typedefs Friends