Rythmos - Transient Integration for Differential Equations Version of the Day
Rythmos_TimeDiscretizedBackwardEulerModelEvaluator.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 
00030 #ifndef RYTHMOS_TIME_DISCRETIZED_BACKWARD_EULER_MODEL_EVALUATOR_HPP
00031 #define RYTHMOS_TIME_DISCRETIZED_BACKWARD_EULER_MODEL_EVALUATOR_HPP
00032 
00033 
00034 #include "Rythmos_Types.hpp"
00035 #include "Thyra_StateFuncModelEvaluatorBase.hpp"
00036 #include "Thyra_ProductVectorBase.hpp"
00037 #include "Thyra_DefaultProductVectorSpace.hpp"
00038 #include "Thyra_DefaultBlockedLinearOp.hpp"
00039 #include "Thyra_DefaultBlockedTriangularLinearOpWithSolveFactory.hpp" // Default implementation
00040 
00041 
00042 namespace Rythmos {
00043 
00044 
00049 template<class Scalar>
00050 class TimeDiscretizedBackwardEulerModelEvaluator
00051   : virtual public Thyra::StateFuncModelEvaluatorBase<Scalar>
00052 {
00053 public:
00054 
00057 
00059   TimeDiscretizedBackwardEulerModelEvaluator();
00060 
00062 
00065 
00067   void initialize(
00068     const RCP<const Thyra::ModelEvaluator<Scalar> > &daeModel,
00069     const Thyra::ModelEvaluatorBase::InArgs<Scalar> &initCond,
00070     const Scalar finalTime,
00071     const int numTimeSteps,
00072     const RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> > &W_bar_factory = Teuchos::null
00073     );
00074 
00076 
00079 
00081   RCP<const Thyra::VectorSpaceBase<Scalar> > get_x_space() const;
00083   RCP<const Thyra::VectorSpaceBase<Scalar> > get_f_space() const;
00085   RCP<Thyra::LinearOpBase<Scalar> > create_W_op() const;
00087   RCP<const Thyra::LinearOpWithSolveFactoryBase<Scalar> > get_W_factory() const;
00089   Thyra::ModelEvaluatorBase::InArgs<Scalar> getNominalValues() const;
00091   Thyra::ModelEvaluatorBase::InArgs<Scalar> createInArgs() const;
00092 
00094 
00095 private:
00096 
00099 
00101   Thyra::ModelEvaluatorBase::OutArgs<Scalar> createOutArgsImpl() const;
00103   void evalModelImpl(
00104     const Thyra::ModelEvaluatorBase::InArgs<Scalar>& inArgs,
00105     const Thyra::ModelEvaluatorBase::OutArgs<Scalar>& outArgs
00106     ) const;
00107 
00109 
00110 private:
00111 
00112   RCP<const Thyra::ModelEvaluator<Scalar> > daeModel_;
00113   Thyra::ModelEvaluatorBase::InArgs<Scalar> initCond_;
00114   Scalar finalTime_;
00115   int numTimeSteps_;
00116 
00117   Scalar initTime_;
00118   Scalar delta_t_;
00119 
00120   RCP<const Thyra::ProductVectorSpaceBase<Scalar> > x_bar_space_;
00121   RCP<const Thyra::ProductVectorSpaceBase<Scalar> > f_bar_space_;
00122   RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> > W_bar_factory_;
00123 
00124 };
00125 
00126 
00131 template<class Scalar>
00132 RCP<TimeDiscretizedBackwardEulerModelEvaluator<Scalar> >
00133 timeDiscretizedBackwardEulerModelEvaluator(
00134   const RCP<const Thyra::ModelEvaluator<Scalar> > &daeModel,
00135   const Thyra::ModelEvaluatorBase::InArgs<Scalar> &initCond,
00136   const Scalar finalTime,
00137   const int numTimeSteps,
00138   const RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> > &W_bar_factory
00139   )
00140 {
00141   RCP<TimeDiscretizedBackwardEulerModelEvaluator<Scalar> >
00142     model(new TimeDiscretizedBackwardEulerModelEvaluator<Scalar>());
00143   model->initialize(daeModel,initCond,finalTime,numTimeSteps,W_bar_factory);
00144   return model;
00145 }
00146 
00147 
00148 // ///////////////////////
00149 // Definition
00150 
00151 
00152 // Constructors/initializers/accessors
00153 
00154 
00155 template<class Scalar>
00156 TimeDiscretizedBackwardEulerModelEvaluator<Scalar>::TimeDiscretizedBackwardEulerModelEvaluator()
00157   :finalTime_(-1.0),
00158    numTimeSteps_(-1),
00159    initTime_(0.0),
00160    delta_t_(-1.0) // Flag for uninitailized!
00161 {}
00162 
00163 
00164 // Overridden from TimeDiscretizedBackwardEulerModelEvaluatorBase
00165 
00166 
00167 template<class Scalar>
00168 void TimeDiscretizedBackwardEulerModelEvaluator<Scalar>::initialize(
00169   const RCP<const Thyra::ModelEvaluator<Scalar> > &daeModel,
00170   const Thyra::ModelEvaluatorBase::InArgs<Scalar> &initCond,
00171   const Scalar finalTime,
00172   const int numTimeSteps,
00173   const RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> > &W_bar_factory
00174   )
00175 {
00176 
00177   TEST_FOR_EXCEPT(is_null(daeModel));
00178   TEST_FOR_EXCEPT(is_null(initCond.get_x()));
00179   TEST_FOR_EXCEPT(is_null(initCond.get_x_dot()));
00180   TEST_FOR_EXCEPT(finalTime <= initCond.get_t());
00181   TEST_FOR_EXCEPT(numTimeSteps <= 0);
00182   // ToDo: Validate that daeModel is of the right form!
00183 
00184   daeModel_ = daeModel;
00185   initCond_ = initCond;
00186   finalTime_ = finalTime;
00187   numTimeSteps_ = numTimeSteps;
00188 
00189   initTime_ = initCond.get_t();
00190   delta_t_ = (finalTime_ - initTime_) / numTimeSteps_;
00191 
00192   x_bar_space_ = productVectorSpace(daeModel_->get_x_space(),numTimeSteps_);
00193   f_bar_space_ = productVectorSpace(daeModel_->get_f_space(),numTimeSteps_);
00194 
00195   if (!is_null(W_bar_factory)) {
00196     W_bar_factory_ = W_bar_factory;
00197   }
00198   else {
00199     W_bar_factory_ =
00200       Thyra::defaultBlockedTriangularLinearOpWithSolveFactory<Scalar>(
00201         daeModel_->get_W_factory()
00202         );
00203   }
00204   
00205 }
00206 
00207 
00208 // Public functions overridden from ModelEvaluator
00209 
00210 
00211 template<class Scalar>
00212 RCP<const Thyra::VectorSpaceBase<Scalar> >
00213 TimeDiscretizedBackwardEulerModelEvaluator<Scalar>::get_x_space() const
00214 {
00215   return x_bar_space_;
00216 }
00217 
00218 
00219 template<class Scalar>
00220 RCP<const Thyra::VectorSpaceBase<Scalar> >
00221 TimeDiscretizedBackwardEulerModelEvaluator<Scalar>::get_f_space() const
00222 {
00223   return f_bar_space_;
00224 }
00225 
00226 
00227 template<class Scalar>
00228 RCP<Thyra::LinearOpBase<Scalar> >
00229 TimeDiscretizedBackwardEulerModelEvaluator<Scalar>::create_W_op() const
00230 {
00231   // Create the block structure for W_op_bar right away!
00232   RCP<Thyra::PhysicallyBlockedLinearOpBase<Scalar> >
00233     W_op_bar = Thyra::defaultBlockedLinearOp<Scalar>();
00234   W_op_bar->beginBlockFill( f_bar_space_, x_bar_space_ );
00235   for ( int k = 0; k < numTimeSteps_; ++k ) {
00236     W_op_bar->setNonconstBlock( k, k, daeModel_->create_W_op() );
00237     if (k > 0)
00238       W_op_bar->setNonconstBlock( k, k-1, daeModel_->create_W_op() );
00239   }
00240   W_op_bar->endBlockFill();
00241   return W_op_bar;
00242 }
00243 
00244 
00245 template<class Scalar>
00246 RCP<const Thyra::LinearOpWithSolveFactoryBase<Scalar> >
00247 TimeDiscretizedBackwardEulerModelEvaluator<Scalar>::get_W_factory() const
00248 {
00249   return W_bar_factory_;
00250 }
00251 
00252 
00253 template<class Scalar>
00254 Thyra::ModelEvaluatorBase::InArgs<Scalar>
00255 TimeDiscretizedBackwardEulerModelEvaluator<Scalar>::getNominalValues() const
00256 {
00257   typedef Thyra::ModelEvaluatorBase MEB;
00258   TEST_FOR_EXCEPT(true);
00259   return MEB::InArgs<Scalar>();
00260 }
00261 
00262 
00263 template<class Scalar>
00264 Thyra::ModelEvaluatorBase::InArgs<Scalar>
00265 TimeDiscretizedBackwardEulerModelEvaluator<Scalar>::createInArgs() const
00266 {
00267   typedef Thyra::ModelEvaluatorBase MEB;
00268   MEB::InArgsSetup<Scalar> inArgs;
00269   inArgs.setModelEvalDescription(this->description());
00270   inArgs.setSupports(MEB::IN_ARG_x);
00271   return inArgs;
00272 }
00273 
00274 
00275 // Private functions overridden from ModelEvaluatorDefaultBase
00276 
00277 
00278 template<class Scalar>
00279 Thyra::ModelEvaluatorBase::OutArgs<Scalar>
00280 TimeDiscretizedBackwardEulerModelEvaluator<Scalar>::createOutArgsImpl() const
00281 {
00282   typedef Thyra::ModelEvaluatorBase MEB;
00283   MEB::OutArgs<Scalar> daeOutArgs = daeModel_->createOutArgs();
00284   MEB::OutArgsSetup<Scalar> outArgs;
00285   outArgs.setModelEvalDescription(this->description());
00286   outArgs.setSupports(MEB::OUT_ARG_f);
00287   outArgs.setSupports(MEB::OUT_ARG_W_op);
00288   outArgs.set_W_properties(daeOutArgs.get_W_properties());
00289   return outArgs;
00290 }
00291 
00292 
00293 template<class Scalar>
00294 void TimeDiscretizedBackwardEulerModelEvaluator<Scalar>::evalModelImpl(
00295   const Thyra::ModelEvaluatorBase::InArgs<Scalar>& inArgs_bar,
00296   const Thyra::ModelEvaluatorBase::OutArgs<Scalar>& outArgs_bar
00297   ) const
00298 {
00299 
00300 
00301   using Teuchos::rcp_dynamic_cast;
00302   typedef ScalarTraits<Scalar> ST;
00303   typedef Thyra::ModelEvaluatorBase MEB;
00304   typedef Thyra::VectorBase<Scalar> VB;
00305   typedef Thyra::ProductVectorBase<Scalar> PVB;
00306   typedef Thyra::BlockedLinearOpBase<Scalar> BLWB;
00307 
00308 /*
00309   THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_GEN_BEGIN(
00310     "Rythmos::ImplicitRKModelEvaluator",inArgs_bar,outArgs_bar,daeModel_
00311     );
00312 */
00313 
00314   TEST_FOR_EXCEPTION( delta_t_ <= 0.0, std::logic_error,
00315     "Error, you have not initialized this object correctly!" );
00316 
00317   //
00318   // A) Unwrap the inArgs and outArgs to get at product vectors and block op
00319   //
00320 
00321   const RCP<const PVB> x_bar = rcp_dynamic_cast<const PVB>(inArgs_bar.get_x(), true);
00322   const RCP<PVB> f_bar = rcp_dynamic_cast<PVB>(outArgs_bar.get_f(), true);
00323   RCP<BLWB> W_op_bar = rcp_dynamic_cast<BLWB>(outArgs_bar.get_W_op(), true);
00324 
00325   //
00326   // B) Assemble f_bar and W_op_bar by looping over stages
00327   //
00328 
00329   MEB::InArgs<Scalar> daeInArgs = daeModel_->createInArgs();
00330   MEB::OutArgs<Scalar> daeOutArgs = daeModel_->createOutArgs();
00331   const RCP<VB> x_dot_i = createMember(daeModel_->get_x_space());
00332   daeInArgs.setArgs(initCond_);
00333   
00334   Scalar t_i = initTime_; // ToDo: Define t_init!
00335 
00336   const Scalar oneOverDeltaT = 1.0/delta_t_;
00337 
00338   for ( int i = 0; i < numTimeSteps_; ++i ) {
00339 
00340     // B.1) Setup the DAE's inArgs for time step eqn f(i) ...
00341     const RCP<const Thyra::VectorBase<Scalar> >
00342       x_i = x_bar->getVectorBlock(i),
00343       x_im1 = ( i==0 ? initCond_.get_x() : x_bar->getVectorBlock(i-1) );
00344     V_VmV( x_dot_i.ptr(), *x_i, *x_im1 ); // x_dot_i = 1/dt * ( x[i] - x[i-1] )
00345     Vt_S( x_dot_i.ptr(), oneOverDeltaT ); // ... 
00346     daeInArgs.set_x_dot( x_dot_i );
00347     daeInArgs.set_x( x_i );
00348     daeInArgs.set_t( t_i );
00349     daeInArgs.set_alpha( oneOverDeltaT );
00350     daeInArgs.set_beta( 1.0 );
00351 
00352     // B.2) Setup the DAE's outArgs for f(i) and/or W(i,i) ...
00353     if (!is_null(f_bar))
00354       daeOutArgs.set_f( f_bar->getNonconstVectorBlock(i) );
00355     if (!is_null(W_op_bar))
00356       daeOutArgs.set_W_op(W_op_bar->getNonconstBlock(i,i).assert_not_null());
00357 
00358     // B.3) Compute f_bar(i) and/or W_op_bar(i,i) ...
00359     daeModel_->evalModel( daeInArgs, daeOutArgs );
00360     daeOutArgs.set_f(Teuchos::null);
00361     daeOutArgs.set_W_op(Teuchos::null);
00362     
00363     // B.4) Evaluate W_op_bar(i,i-1)
00364     if ( !is_null(W_op_bar) && i > 0 ) {
00365       daeInArgs.set_alpha( -oneOverDeltaT );
00366       daeInArgs.set_beta( 0.0 );
00367       daeOutArgs.set_W_op(W_op_bar->getNonconstBlock(i,i-1).assert_not_null());
00368       daeModel_->evalModel( daeInArgs, daeOutArgs );
00369       daeOutArgs.set_W_op(Teuchos::null);
00370     }
00371 
00372     //
00373     t_i += delta_t_;
00374 
00375   }
00376 
00377 /*  
00378   THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_END();
00379 */
00380 
00381 }
00382 
00383 
00384 } // namespace Rythmos
00385 
00386 
00387 #endif // RYTHMOS_TIME_DISCRETIZED_BACKWARD_EULER_MODEL_EVALUATOR_HPP
 All Classes Functions Variables Typedefs Friends