Thyra_DefaultStateEliminationModelEvaluator.hpp

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

Generated on Thu Sep 18 12:33:02 2008 for Thyra Package Browser (Single Doxygen Collection) by doxygen 1.3.9.1