Thyra_DefaultStateEliminationModelEvaluator.hpp

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_DEFAULT_STATE_ELIMINATION_MODEL_EVALUATOR_HPP
00030 #define THYRA_DEFAULT_STATE_ELIMINATION_MODEL_EVALUATOR_HPP
00031 
00032 #include "Thyra_ModelEvaluatorDelegatorBase.hpp"
00033 #include "Thyra_DefaultNominalBoundsOverrideModelEvaluator.hpp"
00034 #include "Thyra_NonlinearSolverBase.hpp"
00035 #include "Teuchos_Time.hpp"
00036 
00037 
00038 namespace Thyra {
00039 
00040 
00048 template<class Scalar>
00049 class DefaultStateEliminationModelEvaluator
00050   : virtual public ModelEvaluatorDelegatorBase<Scalar>
00051 {
00052 public:
00053 
00056 
00058   DefaultStateEliminationModelEvaluator();
00059 
00061   DefaultStateEliminationModelEvaluator(
00062     const Teuchos::RCP<ModelEvaluator<Scalar> > &thyraModel,
00063     const Teuchos::RCP<NonlinearSolverBase<Scalar> > &stateSolver
00064     );
00065 
00067   void initialize(
00068     const Teuchos::RCP<ModelEvaluator<Scalar> > &thyraModel,
00069     const Teuchos::RCP<NonlinearSolverBase<Scalar> > &stateSolver
00070     );
00071 
00073   void uninitialize(
00074     Teuchos::RCP<ModelEvaluator<Scalar> > *thyraModel = NULL,
00075     Teuchos::RCP<NonlinearSolverBase<Scalar> > *stateSolver = NULL
00076     );
00077 
00079 
00082 
00084   std::string description() const;
00085 
00087 
00090 
00092   Teuchos::RCP<const VectorSpaceBase<Scalar> > get_x_space() const;
00094   Teuchos::RCP<const VectorSpaceBase<Scalar> > get_f_space() const;
00096   ModelEvaluatorBase::InArgs<Scalar> getNominalValues() const;
00098   ModelEvaluatorBase::InArgs<Scalar> getLowerBounds() const;
00100   ModelEvaluatorBase::InArgs<Scalar> getUpperBounds() const;
00102   Teuchos::RCP<LinearOpWithSolveBase<Scalar> > create_W() const;
00104   Teuchos::RCP<LinearOpBase<Scalar> > create_W_op() const;
00106   ModelEvaluatorBase::InArgs<Scalar> createInArgs() const;
00107 
00109 
00110 private:
00111 
00114 
00116   ModelEvaluatorBase::OutArgs<Scalar> createOutArgsImpl() const;
00118   void evalModelImpl(
00119     const ModelEvaluatorBase::InArgs<Scalar> &inArgs,
00120     const ModelEvaluatorBase::OutArgs<Scalar> &outArgs
00121     ) const;
00122 
00124 
00125 
00126 private:
00127 
00128   Teuchos::RCP<ModelEvaluator<Scalar> > thyraModel_;
00129   Teuchos::RCP<NonlinearSolverBase<Scalar> > stateSolver_;
00130 
00131   Teuchos::RCP<DefaultNominalBoundsOverrideModelEvaluator<Scalar> > wrappedThyraModel_;
00132 
00133   mutable Teuchos::RCP<VectorBase<Scalar> > x_guess_solu_;
00134   
00135 };
00136 
00137 // /////////////////////////////////
00138 // Implementations
00139 
00140 // Constructors/initializers/accessors/utilities
00141 
00142 template<class Scalar>
00143 DefaultStateEliminationModelEvaluator<Scalar>::DefaultStateEliminationModelEvaluator()
00144 {}
00145 
00146 template<class Scalar>
00147 DefaultStateEliminationModelEvaluator<Scalar>::DefaultStateEliminationModelEvaluator(
00148   const Teuchos::RCP<ModelEvaluator<Scalar> >                 &thyraModel
00149   ,const Teuchos::RCP<NonlinearSolverBase<Scalar> >           &stateSolver
00150   )
00151 {
00152   initialize(thyraModel,stateSolver);
00153 }
00154 
00155 template<class Scalar>
00156 void DefaultStateEliminationModelEvaluator<Scalar>::initialize(
00157   const Teuchos::RCP<ModelEvaluator<Scalar> >                 &thyraModel
00158   ,const Teuchos::RCP<NonlinearSolverBase<Scalar> >           &stateSolver
00159   )
00160 {
00161   this->ModelEvaluatorDelegatorBase<Scalar>::initialize(thyraModel);
00162   TEST_FOR_EXCEPT(!stateSolver.get());
00163   stateSolver_ = stateSolver;
00164   x_guess_solu_ = Teuchos::null; // We will get the guess at the last possible moment!
00165   wrappedThyraModel_ = Teuchos::rcp(
00166     new DefaultNominalBoundsOverrideModelEvaluator<Scalar>(
00167       Teuchos::rcp_const_cast<ModelEvaluator<Scalar> >(thyraModel)
00168       ,Teuchos::null
00169       )
00170     );
00171   stateSolver_->setModel(wrappedThyraModel_);
00172 }
00173 
00174 template<class Scalar>
00175 void DefaultStateEliminationModelEvaluator<Scalar>::uninitialize(
00176   Teuchos::RCP<ModelEvaluator<Scalar> >                 *thyraModel
00177   ,Teuchos::RCP<NonlinearSolverBase<Scalar> >           *stateSolver
00178   )
00179 {
00180   if(thyraModel) *thyraModel = this->getUnderlyingModel();
00181   if(stateSolver) *stateSolver = stateSolver_;
00182   this->ModelEvaluatorDelegatorBase<Scalar>::uninitialize();
00183   stateSolver_ = Teuchos::null;
00184   wrappedThyraModel_ = Teuchos::null;
00185   x_guess_solu_ = Teuchos::null;
00186 }
00187 
00188 // Public functions overridden from Teuchos::Describable
00189 
00190 
00191 template<class Scalar>
00192 std::string DefaultStateEliminationModelEvaluator<Scalar>::description() const
00193 {
00194   const Teuchos::RCP<const ModelEvaluator<Scalar> >
00195     thyraModel = this->getUnderlyingModel();
00196   std::ostringstream oss;
00197   oss << "Thyra::DefaultStateEliminationModelEvaluator{";
00198   oss << "thyraModel=";
00199   if(thyraModel.get())
00200     oss << "\'"<<thyraModel->description()<<"\'";
00201   else
00202     oss << "NULL";
00203   oss << ",stateSolver=";
00204   if(stateSolver_.get())
00205     oss << "\'"<<stateSolver_->description()<<"\'";
00206   else
00207     oss << "NULL";
00208   oss << "}";
00209   return oss.str();
00210 }
00211 
00212 // Public functions overridden from ModelEvaulator
00213 
00214 template<class Scalar>
00215 Teuchos::RCP<const VectorSpaceBase<Scalar> >
00216 DefaultStateEliminationModelEvaluator<Scalar>::get_x_space() const
00217 {
00218   return Teuchos::null;
00219 }
00220 
00221 template<class Scalar>
00222 Teuchos::RCP<const VectorSpaceBase<Scalar> >
00223 DefaultStateEliminationModelEvaluator<Scalar>::get_f_space() const
00224 {
00225   return Teuchos::null;
00226 }
00227 
00228 template<class Scalar>
00229 ModelEvaluatorBase::InArgs<Scalar>
00230 DefaultStateEliminationModelEvaluator<Scalar>::getNominalValues() const
00231 {
00232   typedef ModelEvaluatorBase MEB;
00233   const Teuchos::RCP<const ModelEvaluator<Scalar> >
00234     thyraModel = this->getUnderlyingModel();
00235   MEB::InArgsSetup<Scalar> nominalValues(thyraModel->getNominalValues());
00236   nominalValues.setModelEvalDescription(this->description());
00237   nominalValues.setUnsupportsAndRelated(MEB::IN_ARG_x); // Wipe out x, x_dot ...
00238   return nominalValues;
00239 }
00240 
00241 template<class Scalar>
00242 ModelEvaluatorBase::InArgs<Scalar>
00243 DefaultStateEliminationModelEvaluator<Scalar>::getLowerBounds() const
00244 {
00245   typedef ModelEvaluatorBase MEB;
00246   const Teuchos::RCP<const ModelEvaluator<Scalar> >
00247     thyraModel = this->getUnderlyingModel();
00248   MEB::InArgsSetup<Scalar> lowerBounds(thyraModel->getLowerBounds());
00249   lowerBounds.setModelEvalDescription(this->description());
00250   lowerBounds.setUnsupportsAndRelated(MEB::IN_ARG_x); // Wipe out x, x_dot ...
00251   return lowerBounds;
00252 }
00253 
00254 template<class Scalar>
00255 ModelEvaluatorBase::InArgs<Scalar>
00256 DefaultStateEliminationModelEvaluator<Scalar>::getUpperBounds() const
00257 {
00258   typedef ModelEvaluatorBase MEB;
00259   const Teuchos::RCP<const ModelEvaluator<Scalar> >
00260     thyraModel = this->getUnderlyingModel();
00261   MEB::InArgsSetup<Scalar> upperBounds(thyraModel->getUpperBounds());
00262   upperBounds.setModelEvalDescription(this->description());
00263   upperBounds.setUnsupportsAndRelated(MEB::IN_ARG_x); // Wipe out x, x_dot ...
00264   return upperBounds;
00265 }
00266 
00267 template<class Scalar>
00268 Teuchos::RCP<LinearOpWithSolveBase<Scalar> >
00269 DefaultStateEliminationModelEvaluator<Scalar>::create_W() const
00270 {
00271   return Teuchos::null;
00272 }
00273 
00274 template<class Scalar>
00275 Teuchos::RCP<LinearOpBase<Scalar> >
00276 DefaultStateEliminationModelEvaluator<Scalar>::create_W_op() const
00277 {
00278   return Teuchos::null;
00279 }
00280 
00281 
00282 template<class Scalar>
00283 ModelEvaluatorBase::InArgs<Scalar>
00284 DefaultStateEliminationModelEvaluator<Scalar>::createInArgs() const
00285 {
00286   typedef ModelEvaluatorBase MEB;
00287   const Teuchos::RCP<const ModelEvaluator<Scalar> >
00288     thyraModel = this->getUnderlyingModel();
00289   const MEB::InArgs<Scalar> wrappedInArgs = thyraModel->createInArgs();
00290   MEB::InArgsSetup<Scalar> inArgs;
00291   inArgs.setModelEvalDescription(this->description());
00292   inArgs.set_Np(wrappedInArgs.Np());
00293   inArgs.setSupports(wrappedInArgs);
00294   inArgs.setUnsupportsAndRelated(MEB::IN_ARG_x); // Wipe out x, x_dot ...
00295   return inArgs;
00296 }
00297 
00298 
00299 // Private functions overridden from ModelEvaulatorDefaultBase
00300 
00301 
00302 template<class Scalar>
00303 ModelEvaluatorBase::OutArgs<Scalar>
00304 DefaultStateEliminationModelEvaluator<Scalar>::createOutArgsImpl() const
00305 {
00306   typedef ModelEvaluatorBase MEB;
00307   const Teuchos::RCP<const ModelEvaluator<Scalar> >
00308     thyraModel = this->getUnderlyingModel();
00309   const MEB::OutArgs<Scalar> wrappedOutArgs = thyraModel->createOutArgs();
00310   const int Np = wrappedOutArgs.Np(), Ng = wrappedOutArgs.Ng();
00311   MEB::OutArgsSetup<Scalar> outArgs;
00312   outArgs.setModelEvalDescription(this->description());
00313   outArgs.set_Np_Ng(Np,Ng);
00314   outArgs.setSupports(wrappedOutArgs);
00315   outArgs.setUnsupportsAndRelated(MEB::IN_ARG_x); // wipe out DgDx ...
00316   outArgs.setUnsupportsAndRelated(MEB::OUT_ARG_f); // wipe out f, W, DfDp ...
00317   return outArgs;
00318 }
00319 
00320 template<class Scalar>
00321 void DefaultStateEliminationModelEvaluator<Scalar>::evalModelImpl(
00322   const ModelEvaluatorBase::InArgs<Scalar> &inArgs,
00323   const ModelEvaluatorBase::OutArgs<Scalar> &outArgs
00324   ) const
00325 {
00326   typedef ModelEvaluatorBase MEB;
00327   using Teuchos::RCP;
00328   using Teuchos::rcp;
00329   using Teuchos::rcp_const_cast;
00330   using Teuchos::rcp_dynamic_cast;
00331   using Teuchos::OSTab;
00332 
00333   Teuchos::Time totalTimer(""), timer("");
00334   totalTimer.start(true);
00335 
00336   const Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
00337   const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00338   Teuchos::OSTab tab(out);
00339   if(out.get() && static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_LOW))
00340     *out << "\nEntering Thyra::DefaultStateEliminationModelEvaluator<Scalar>::evalModel(...) ...\n";
00341 
00342   const Teuchos::RCP<const ModelEvaluator<Scalar> >
00343     thyraModel = this->getUnderlyingModel();
00344 
00345   const int Np = outArgs.Np(), Ng = outArgs.Ng();
00346 
00347   // Get the intial state guess if not already gotten
00348   if (is_null(x_guess_solu_)) {
00349     const ModelEvaluatorBase::InArgs<Scalar>
00350       nominalValues = thyraModel->getNominalValues();
00351     if(nominalValues.get_x().get()) {
00352       x_guess_solu_ = nominalValues.get_x()->clone_v();
00353     }
00354     else {
00355       x_guess_solu_ = createMember(thyraModel->get_x_space());
00356       assign(&*x_guess_solu_,Scalar(0.0));
00357     }
00358   }
00359 
00360   // Reset the nominal values
00361   MEB::InArgs<Scalar> wrappedNominalValues = thyraModel->getNominalValues();
00362   wrappedNominalValues.setArgs(inArgs,true);
00363   wrappedNominalValues.set_x(x_guess_solu_);
00364   
00365   typedef Teuchos::VerboseObjectTempState<ModelEvaluatorBase> VOTSME;
00366   //VOTSME thyraModel_outputTempState(rcp(&wrappedThyraModel,false),out,verbLevel);
00367 
00368   typedef Teuchos::VerboseObjectTempState<NonlinearSolverBase<Scalar> > VOTSNSB;
00369   VOTSNSB statSolver_outputTempState(
00370     stateSolver_,out
00371     ,static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_LOW) ? Teuchos::VERB_LOW : Teuchos::VERB_NONE 
00372     );
00373 
00374   if(out.get() && static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_EXTREME))
00375     *out
00376       << "\ninArgs =\n" << Teuchos::describe(inArgs,verbLevel)
00377       << "\noutArgs on input =\n" << Teuchos::describe(outArgs,Teuchos::VERB_LOW);
00378 
00379   if(out.get() && static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_LOW))
00380     *out << "\nSolving f(x,...) for x ...\n";
00381 
00382   wrappedThyraModel_->setNominalValues(
00383     rcp(new MEB::InArgs<Scalar>(wrappedNominalValues))
00384     );
00385   
00386   SolveStatus<Scalar> solveStatus = stateSolver_->solve(&*x_guess_solu_,NULL);
00387 
00388   if( solveStatus.solveStatus == SOLVE_STATUS_CONVERGED ) {
00389     
00390     if(out.get() && static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_LOW))
00391       *out << "\nComputing the output functions at the solved state solution ...\n";
00392 
00393     MEB::InArgs<Scalar>   wrappedInArgs  = thyraModel->createInArgs();
00394     MEB::OutArgs<Scalar>  wrappedOutArgs = thyraModel->createOutArgs();
00395     wrappedInArgs.setArgs(inArgs,true);
00396     wrappedInArgs.set_x(x_guess_solu_);
00397     wrappedOutArgs.setArgs(outArgs,true);
00398     
00399     for( int l = 0; l < Np; ++l ) {
00400       for( int j = 0; j < Ng; ++j ) {
00401         if(
00402           outArgs.supports(MEB::OUT_ARG_DgDp,j,l).none()==false
00403           && outArgs.get_DgDp(j,l).isEmpty()==false
00404           )
00405         {
00406           // Set DfDp(l) and DgDx(j) to be computed!
00407           //wrappedOutArgs.set_DfDp(l,...);
00408           //wrappedOutArgs.set_DgDx(j,...);
00409           TEST_FOR_EXCEPT(true);
00410         }
00411       }
00412     }
00413     
00414     thyraModel->evalModel(wrappedInArgs,wrappedOutArgs);
00415 
00416     //
00417     // Compute DgDp(j,l) using direct sensitivties
00418     //
00419     for( int l = 0; l < Np; ++l ) {
00420       if(
00421         wrappedOutArgs.supports(MEB::OUT_ARG_DfDp,l).none()==false
00422         && wrappedOutArgs.get_DfDp(l).isEmpty()==false
00423         )
00424       {
00425         //
00426         // Compute:  D(l) = -inv(DfDx)*DfDp(l)
00427         //
00428         TEST_FOR_EXCEPT(true);
00429         for( int j = 0; j < Ng; ++j ) {
00430           if(
00431             outArgs.supports(MEB::OUT_ARG_DgDp,j,l).none()==false
00432             && outArgs.get_DgDp(j,l).isEmpty()==false
00433             )
00434           {
00435             //
00436             // Compute:  DgDp(j,l) = DgDp(j,l) + DgDx(j)*D
00437             //
00438             TEST_FOR_EXCEPT(true);
00439           }
00440         }
00441       }
00442     }
00443     // ToDo: Add a mode to compute DgDp(l) using adjoint sensitivities?
00444     
00445   }
00446   else {
00447     
00448     if(out.get() && static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_LOW))
00449       *out << "\nFailed to converge, returning NaNs ...\n";
00450     outArgs.setFailed();
00451     
00452   }
00453   
00454   if(out.get() && static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_EXTREME))
00455     *out
00456       << "\noutArgs on output =\n" << Teuchos::describe(outArgs,verbLevel);
00457 
00458   totalTimer.stop();
00459   if(out.get() && static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_LOW))
00460     *out
00461       << "\nTotal evaluation time = "<<totalTimer.totalElapsedTime()<<" sec\n"
00462       << "\nLeaving Thyra::DefaultStateEliminationModelEvaluator<Scalar>::evalModel(...) ...\n";
00463   
00464 }
00465 
00466 
00467 } // namespace Thyra
00468 
00469 
00470 #endif // THYRA_DEFAULT_STATE_ELIMINATION_MODEL_EVALUATOR_HPP

Generated on Tue Jul 13 09:26:21 2010 for Thyra by  doxygen 1.4.7