Thyra Package Browser (Single Doxygen Collection): Thyra_SingleResidSSDAEModelEvaluator.hpp Source File

00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 //    Thyra: Interfaces and Support for Abstract Numerical Algorithms
00005 //                 Copyright (2004) 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 Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #ifndef THYRA_SINGLE_RESID_SS_DAE_MODEL_EVALUATOR_HPP
00030 #define THYRA_SINGLE_RESID_SS_DAE_MODEL_EVALUATOR_HPP
00031 
00032 #include "Thyra_ModelEvaluatorDelegatorBase.hpp"
00033 #include "Thyra_ModelEvaluatorHelpers.hpp"
00034 #include "Thyra_VectorStdOps.hpp"
00035 
00036 #ifdef THYRA_RYTHMOS_DEBUG
00037 #include "Thyra_TestingTools.hpp"
00038 #endif // THYRA_RYTHMOS_DEBUG
00039 
00040 namespace Thyra {
00041 
00063 template<class Scalar>
00064 class SingleResidSSDAEModelEvaluator
00065   : virtual public ModelEvaluatorDelegatorBase<Scalar>
00066 {
00067 public:
00068 
00071 
00073   SingleResidSSDAEModelEvaluator();
00074 
00076   SingleResidSSDAEModelEvaluator(
00077     const Teuchos::RefCountPtr<const ModelEvaluator<Scalar> >     &daeModel
00078     ,const Scalar                                                 &coeff_x_dot
00079     ,const Teuchos::RefCountPtr<const VectorBase<Scalar> >        &x_dot_base
00080     ,const Scalar                                                 &coeff_x
00081     ,const Teuchos::RefCountPtr<const VectorBase<Scalar> >        &x_base
00082     ,const Scalar                                                 &t_base
00083     ,const Teuchos::RefCountPtr<const VectorBase<Scalar> >        &x_bar_init
00084     );
00085 
00087   void initialize(
00088     const Teuchos::RefCountPtr<const ModelEvaluator<Scalar> >     &daeModel
00089     ,const Scalar                                                 &coeff_x_dot
00090     ,const Teuchos::RefCountPtr<const VectorBase<Scalar> >        &x_dot_base
00091     ,const Scalar                                                 &coeff_x
00092     ,const Teuchos::RefCountPtr<const VectorBase<Scalar> >        &x_base
00093     ,const Scalar                                                 &t_base
00094     ,const Teuchos::RefCountPtr<const VectorBase<Scalar> >        &x_bar_init
00095     );
00096 
00098 
00101 
00103   ModelEvaluatorBase::InArgs<Scalar> getNominalValues() const;
00105   ModelEvaluatorBase::InArgs<Scalar> createInArgs() const;
00107   ModelEvaluatorBase::OutArgs<Scalar> createOutArgs() const;
00109   void evalModel( const ModelEvaluatorBase::InArgs<Scalar>& inArgs, const ModelEvaluatorBase::OutArgs<Scalar>& outArgs ) const;
00110 
00112 
00113 private:
00114 
00115   Scalar                                                 coeff_x_dot_;
00116   Teuchos::RefCountPtr<const VectorBase<Scalar> >        x_dot_base_;
00117   Scalar                                                 coeff_x_;
00118   Teuchos::RefCountPtr<const VectorBase<Scalar> >        x_base_;
00119   Scalar                                                 t_base_;
00120 
00121   ModelEvaluatorBase::InArgs<Scalar>                     nominalValues_;
00122 
00123   // cache
00124   Teuchos::RefCountPtr<VectorBase<Scalar> >        x_;
00125   Teuchos::RefCountPtr<VectorBase<Scalar> >        x_dot_;
00126 
00127 
00128 };
00129 
00130 // ///////////////////////
00131 // Definition
00132 
00133 // Constructors/initializers/accessors
00134 
00135 template<class Scalar>
00136 SingleResidSSDAEModelEvaluator<Scalar>::SingleResidSSDAEModelEvaluator()
00137 {}
00138 
00139 template<class Scalar>
00140 SingleResidSSDAEModelEvaluator<Scalar>::SingleResidSSDAEModelEvaluator(
00141   const Teuchos::RefCountPtr<const ModelEvaluator<Scalar> >     &daeModel
00142   ,const Scalar                                                 &coeff_x_dot
00143   ,const Teuchos::RefCountPtr<const VectorBase<Scalar> >        &x_dot_base
00144   ,const Scalar                                                 &coeff_x
00145   ,const Teuchos::RefCountPtr<const VectorBase<Scalar> >        &x_base
00146   ,const Scalar                                                 &t_base
00147   ,const Teuchos::RefCountPtr<const VectorBase<Scalar> >        &x_bar_init
00148   )
00149 {
00150   initialize(daeModel,coeff_x_dot,x_dot_base,coeff_x,x_base,t_base,x_bar_init);
00151 }
00152 
00153 template<class Scalar>
00154 void SingleResidSSDAEModelEvaluator<Scalar>::initialize(
00155   const Teuchos::RefCountPtr<const ModelEvaluator<Scalar> >     &daeModel
00156   ,const Scalar                                                 &coeff_x_dot
00157   ,const Teuchos::RefCountPtr<const VectorBase<Scalar> >        &x_dot_base
00158   ,const Scalar                                                 &coeff_x
00159   ,const Teuchos::RefCountPtr<const VectorBase<Scalar> >        &x_base
00160   ,const Scalar                                                 &t_base
00161   ,const Teuchos::RefCountPtr<const VectorBase<Scalar> >        &x_bar_init
00162   )
00163 {
00164   this->ModelEvaluatorDelegatorBase<Scalar>::initialize(daeModel);
00165   coeff_x_dot_   = coeff_x_dot;
00166   x_dot_base_    = x_dot_base;
00167   coeff_x_       = coeff_x;
00168   x_base_        = x_base;
00169   t_base_        = t_base;
00170 
00171   nominalValues_ = daeModel->getNominalValues();
00172   nominalValues_.set_x(x_bar_init);
00173 
00174   x_dot_ = createMember( daeModel->get_x_space() );
00175   x_ = createMember( daeModel->get_x_space() );
00176 
00177   // ToDo: Check that daeModel supports x_dot, x and maybe t
00178 
00179 #ifdef THYRA_RYTHMOS_DEBUG
00180   std::cout << "----------------------------------------------------------------------" << std::endl;
00181   std::cout << "Thyra::SingleResidSSDAEModelEvaluator::initialize" << std::endl;
00182   std::cout << "coeff_x_dot_ = " << coeff_x_dot_ << std::endl;
00183   std::cout << "x_dot_base_ = ";                  
00184   if ( x_dot_base_.get() ) 
00185     std::cout << "\n" << *x_dot_base_            << std::endl;
00186   else
00187     std::cout << "null"                          << std::endl;
00188   std::cout << "coeff_x_ = " << coeff_x_         << std::endl;
00189   std::cout << "x_base_ = ";                      
00190   if ( x_base_.get() )
00191     std::cout << "\n" << *x_base_                << std::endl;
00192   else
00193     std::cout << "null"                          << std::endl;
00194   std::cout << "t_base_ = " << t_base_           << std::endl;
00195   std::cout << "x_bar_init = ";                  
00196   if ( x_bar_init.get() )
00197     std::cout << "\n" <<  *x_bar_init_           << std::endl;
00198   else
00199     std::cout << "null"                          << std::endl;
00200   std::cout << "x_dot_ = ";                       
00201   if ( x_dot_.get() )
00202     std::cout << "\n" << *x_dot_                 << std::endl;
00203   else
00204     std::cout << "null"                          << std::endl;
00205   std::cout << "x_ = ";                           
00206   if ( x_.get() )
00207     std::cout << "\n" << *x_                     << std::endl;
00208   else
00209     std::cout << "null"                          << std::endl;
00210   std::cout << "----------------------------------------------------------------------" << std::endl;
00211 #endif // THYRA_RYTHMOS_DEBUG
00212 }
00213 
00214 // Overridden from ModelEvaluator
00215 
00216 template<class Scalar>
00217 ModelEvaluatorBase::InArgs<Scalar>
00218 SingleResidSSDAEModelEvaluator<Scalar>::getNominalValues() const
00219 {
00220   return nominalValues_;
00221 }
00222 
00223 template<class Scalar>
00224 ModelEvaluatorBase::InArgs<Scalar>
00225 SingleResidSSDAEModelEvaluator<Scalar>::createInArgs() const
00226 {
00227   typedef ModelEvaluatorBase MEB;
00228   MEB::InArgsSetup<Scalar> inArgs;
00229   inArgs.setSupports(MEB::IN_ARG_x);
00230   inArgs.setSupports(MEB::IN_ARG_beta);
00231   return inArgs;
00232 }
00233 
00234 template<class Scalar>
00235 ModelEvaluatorBase::OutArgs<Scalar>
00236 SingleResidSSDAEModelEvaluator<Scalar>::createOutArgs() const
00237 {
00238   typedef ModelEvaluatorBase MEB;
00239   MEB::OutArgsSetup<Scalar> outArgs;
00240   outArgs.setSupports(MEB::OUT_ARG_f);
00241   outArgs.setSupports(MEB::OUT_ARG_W);
00242   return outArgs;
00243 }
00244 
00245 template<class Scalar>
00246 void SingleResidSSDAEModelEvaluator<Scalar>::evalModel(
00247   const ModelEvaluatorBase::InArgs<Scalar>&   inArgs_bar
00248   ,const ModelEvaluatorBase::OutArgs<Scalar>& outArgs_bar
00249   ) const
00250 {
00251   const Teuchos::RefCountPtr<const ModelEvaluator<Scalar> >
00252     daeModel = this->getUnderlyingModel();
00253 #ifdef THYRA_RYTHMOS_DEBUG
00254   std::cout << "----------------------------------------------------------------------" << std::endl;
00255   std::cout << "Thyra::SingleResidSSDAEModelEvaluator::evalModel" << std::endl;
00256 #endif // THYRA_RYTHMOS_DEBUG
00257   const VectorBase<Scalar> &x_bar = *inArgs_bar.get_x(); 
00258   // x_dot = coeff_x_dot * x_bar + x_dot_base
00259   if (x_dot_base_.get())
00260     Thyra::V_StVpV( &*x_dot_, coeff_x_dot_, x_bar, *x_dot_base_ );
00261   else
00262     Thyra::V_StV( &*x_dot_, coeff_x_dot_, x_bar);
00263 #ifdef THYRA_RYTHMOS_DEBUG
00264   std::cout << "x_dot_ = coeff_x_dot_ * x_bar + x_dot_base_" << std::endl;
00265   std::cout << "coeff_x_dot_ = " << coeff_x_dot_             << std::endl;
00266   std::cout << "x_bar = "                                    << std::endl;
00267   std::cout <<  x_bar                                        << std::endl;
00268   std::cout << "x_dot_base_ = ";
00269   if ( x_dot_base_.get() )
00270     std::cout << "\n" << *x_dot_base_                        << std::endl;
00271   else
00272     std::cout << "null"                                      << std::endl;
00273   std::cout << "x_dot_ = ";
00274   if ( x_dot_.get() )
00275     std::cout << "\n" << *x_dot_                             << std::endl;
00276   else
00277     std::cout << "null"                                      << std::endl;
00278 #endif // THYRA_RYTHMOS_DEBUG
00279 
00280   // x = coeff_x * x_bar + x_base
00281   if (x_base_.get())
00282     Thyra::V_StVpV( &*x_, coeff_x_, x_bar, *x_base_ );
00283   else
00284     Thyra::V_StV( &*x_, coeff_x_, x_bar);
00285 #ifdef THYRA_RYTHMOS_DEBUG
00286   std::cout << "x_ = coeff_x_ * x_bar + x_base_" << std::endl;
00287   std::cout << "coeff_x_ = " << coeff_x_         << std::endl;
00288   std::cout << "x_bar = "                        << std::endl;
00289   std::cout <<  x_bar                            << std::endl;
00290   std::cout << "x_base_ = ";
00291   if ( x_base_.get() )
00292     std::cout << "\n" << *x_base_                << std::endl;
00293   else
00294     std::cout << "null"                          << std::endl;
00295   std::cout << "x_ = ";
00296   if ( x_.get() )
00297     std::cout << "\n" << *x_                     << std::endl;
00298   else
00299     std::cout << "null"                          << std::endl;
00300 #endif // THYRA_RYTHMOS_DEBUG
00301 
00302   // Compute W and f
00303   Teuchos::RefCountPtr<LinearOpWithSolveBase<Scalar> > W;
00304   if( (W = outArgs_bar.get_W()).get() ) {
00305     // Compute Jacobian and the residual
00306     const Scalar beta = inArgs_bar.get_beta();
00307     eval_f_W(
00308       *daeModel
00309       ,*x_dot_, *x_, t_base_, Scalar(beta*coeff_x_dot_), Scalar(beta*coeff_x_)
00310       ,outArgs_bar.get_f().get(), &*W
00311       );
00312 #ifdef THYRA_RYTHMOS_DEBUG
00313     std::cout << "f = "                 << std::endl;
00314     std::cout << *(outArgs_bar.get_f()) << std::endl;
00315     std::cout << "W = "                 << std::endl;
00316     std::cout << *W                     << std::endl;
00317 #endif // THYRA_RYTHMOS_DEBUG
00318   }
00319   else {
00320     // Compute only the residual
00321     eval_f( *daeModel, *x_dot_, *x_, t_base_, &*outArgs_bar.get_f() );
00322 #ifdef THYRA_RYTHMOS_DEBUG
00323     std::cout << "f = "                 << std::endl;
00324     std::cout << *(outArgs_bar.get_f()) << std::endl;
00325 #endif // THYRA_RYTHMOS_DEBUG
00326   }
00327 #ifdef THYRA_RYTHMOS_DEBUG
00328   std::cout << "----------------------------------------------------------------------" << std::endl;
00329 #endif // THYRA_RYTHMOS_DEBUG
00330 }
00331 
00332 } // namespace Thyra
00333 
00334 #endif // THYRA_SINGLE_RESID_SS_DAE_MODEL_EVALUATOR_HPP