EpetraExt Package Browser (Single Doxygen Collection) Development
GLpApp_AdvDiffReactOptModel.cpp
Go to the documentation of this file.
00001 /*
00002 //@HEADER
00003 // ***********************************************************************
00004 //
00005 //     EpetraExt: Epetra Extended - Linear Algebra Services Package
00006 //                 Copyright (2011) Sandia Corporation
00007 //
00008 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00009 // the U.S. Government retains certain rights in this software.
00010 //
00011 // Redistribution and use in source and binary forms, with or without
00012 // modification, are permitted provided that the following conditions are
00013 // met:
00014 //
00015 // 1. Redistributions of source code must retain the above copyright
00016 // notice, this list of conditions and the following disclaimer.
00017 //
00018 // 2. Redistributions in binary form must reproduce the above copyright
00019 // notice, this list of conditions and the following disclaimer in the
00020 // documentation and/or other materials provided with the distribution.
00021 //
00022 // 3. Neither the name of the Corporation nor the names of the
00023 // contributors may be used to endorse or promote products derived from
00024 // this software without specific prior written permission.
00025 //
00026 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00027 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00028 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00029 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00030 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00031 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00032 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00033 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00034 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00035 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00036 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00037 //
00038 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00039 //
00040 // ***********************************************************************
00041 //@HEADER
00042 */
00043 
00044 #include "GLpApp_AdvDiffReactOptModel.hpp"
00045 #include "Epetra_SerialComm.h"
00046 #include "Epetra_CrsMatrix.h"
00047 #include "Teuchos_ScalarTraits.hpp"
00048 #include "Teuchos_VerboseObject.hpp"
00049 #include "Teuchos_TimeMonitor.hpp"
00050 
00051 // 2007/06/14: rabartl: The dependence on Thyra is non-optional right now
00052 // since I use it to perform a mat-vec that I could not figure out how to do
00053 // with just epetra.  If this becomes a problem then we have to change this
00054 // later!  Just grep for 'Thyra' outside of the ifdef for
00055 // HAVE_THYRA_EPETRAEXT.
00056 #include "Thyra_EpetraThyraWrappers.hpp"
00057 #include "Thyra_VectorStdOps.hpp"
00058 
00059 #ifdef HAVE_THYRA_EPETRAEXT
00060 // For orthogonalization of the basis B_bar
00061 #include "sillyModifiedGramSchmidt.hpp" // This is just an example!
00062 #include "Thyra_MultiVectorStdOps.hpp"
00063 #include "Thyra_SpmdVectorSpaceBase.hpp"
00064 #endif // HAVE_THYRA_EPETRAEXT
00065 
00066 //#define GLPAPP_ADVDIFFREACT_OPTMODEL_DUMP_STUFF
00067 
00068 namespace {
00069 
00070 Teuchos::RefCountPtr<Teuchos::Time>
00071   initalizationTimer    = Teuchos::TimeMonitor::getNewTimer("AdvDiffReact:Init"),
00072   evalTimer             = Teuchos::TimeMonitor::getNewTimer("AdvDiffReact:Eval"),
00073   p_bar_Timer           = Teuchos::TimeMonitor::getNewTimer("AdvDiffReact:Eval:p_bar"),
00074   R_p_bar_Timer         = Teuchos::TimeMonitor::getNewTimer("AdvDiffReact:Eval:R_p_bar"),
00075   f_Timer               = Teuchos::TimeMonitor::getNewTimer("AdvDiffReact:Eval:f"),
00076   g_Timer               = Teuchos::TimeMonitor::getNewTimer("AdvDiffReact:Eval:g"),
00077   W_Timer               = Teuchos::TimeMonitor::getNewTimer("AdvDiffReact:Eval:W"),
00078   DfDp_Timer            = Teuchos::TimeMonitor::getNewTimer("AdvDiffReact:Eval:DfDp"),
00079   DgDx_Timer            = Teuchos::TimeMonitor::getNewTimer("AdvDiffReact:Eval:DgDx"),
00080   DgDp_Timer            = Teuchos::TimeMonitor::getNewTimer("AdvDiffReact:Eval:DgDp");
00081 
00082 inline double sqr( const double& s ) { return s*s; }
00083 
00084 inline double dot( const Epetra_Vector &x, const Epetra_Vector &y )
00085 {
00086   double dot[1];
00087   x.Dot(y,dot);
00088   return dot[0];
00089 }
00090 
00091 } // namespace
00092 
00093 namespace GLpApp {
00094 
00095 AdvDiffReactOptModel::AdvDiffReactOptModel(
00096   const Teuchos::RefCountPtr<const Epetra_Comm>  &comm
00097   ,const double                                  beta
00098   ,const double                                  len_x     // Ignored if meshFile is *not* empty
00099   ,const double                                  len_y     // Ignored if meshFile is *not* empty
00100   ,const int                                     local_nx  // Ignored if meshFile is *not* empty
00101   ,const int                                     local_ny  // Ignored if meshFile is *not* empty
00102   ,const char                                    meshFile[]
00103   ,const int                                     np
00104   ,const double                                  x0
00105   ,const double                                  p0
00106   ,const double                                  reactionRate
00107   ,const bool                                    normalizeBasis
00108   ,const bool                                    supportDerivatives
00109   )
00110   : supportDerivatives_(supportDerivatives), np_(np)
00111 {
00112   Teuchos::TimeMonitor initalizationTimerMonitor(*initalizationTimer);
00113 #ifdef GLPAPP_ADVDIFFREACT_OPTMODEL_DUMP_STUFF
00114   Teuchos::RefCountPtr<Teuchos::FancyOStream>
00115     out = Teuchos::VerboseObjectBase::getDefaultOStream();
00116   Teuchos::OSTab tab(out);
00117   *out << "\nEntering AdvDiffReactOptModel::AdvDiffReactOptModel(...) ...\n";
00118 #endif
00119   //
00120   // Initalize the data pool object
00121   //
00122   dat_ = Teuchos::rcp(new GLpApp::GLpYUEpetraDataPool(comm,beta,len_x,len_y,local_nx,local_ny,meshFile,false));
00123   //
00124   // Get the maps
00125   //
00126   map_x_ = Teuchos::rcp(new Epetra_Map(dat_->getA()->OperatorDomainMap()));
00127   map_p_.resize(Np_);
00128   map_p_[p_rx_idx] = Teuchos::rcp(new Epetra_Map(1,1,0,*comm));
00129   map_p_bar_ = Teuchos::rcp(new Epetra_Map(dat_->getB()->OperatorDomainMap()));
00130   map_f_ = Teuchos::rcp(new Epetra_Map(dat_->getA()->OperatorRangeMap()));
00131   map_g_ = Teuchos::rcp(new Epetra_Map(1,1,0,*comm));
00132   //
00133   // Initialize the basis matrix for p such that p_bar = B_bar * p
00134   //
00135   if( np_ > 0 ) {
00136     //
00137     // Create a sine series basis for y (the odd part of a Fourier series basis)
00138     //
00139     const Epetra_SerialDenseMatrix &ipcoords = *dat_->getipcoords();
00140     const Epetra_IntSerialDenseVector &pindx = *dat_->getpindx();
00141     Epetra_Map overlapmap(-1,pindx.M(),const_cast<int*>(pindx.A()),1,*comm);
00142     const double pi = 2.0 * std::asin(1.0);
00143 #ifdef GLPAPP_ADVDIFFREACT_OPTMODEL_DUMP_STUFF
00144     *out << "\npi="<<pi<<"\n";
00145 #endif
00146     map_p_[p_bndy_idx] = Teuchos::rcp(new Epetra_Map(np_,np_,0,*comm));
00147     B_bar_ = Teuchos::rcp(new Epetra_MultiVector(*map_p_bar_,np_));
00148     (*B_bar_)(0)->PutScalar(1.0); // First column is all ones!
00149     if( np_ > 1 ) {
00150       const int numBndyNodes        = map_p_bar_->NumMyElements();
00151       const int *bndyNodeGlobalIDs  = map_p_bar_->MyGlobalElements();
00152       for( int i = 0; i < numBndyNodes; ++i ) {
00153 #ifdef GLPAPP_ADVDIFFREACT_OPTMODEL_DUMP_STUFF
00154         const int global_id = bndyNodeGlobalIDs[i];
00155 #endif
00156         const int local_id = overlapmap.LID(bndyNodeGlobalIDs[i]);
00157         const double x = ipcoords(local_id,0), y = ipcoords(local_id,1);
00158         double z = -1.0, L = -1.0;
00159         if( x < 1e-10 || len_x - 1e-10 < x ) {
00160           z = y;
00161           L = len_y;
00162         }
00163         else {
00164           z = x;
00165           L = len_x;
00166         }
00167 #ifdef GLPAPP_ADVDIFFREACT_OPTMODEL_DUMP_STUFF
00168         *out << "\ni="<<i<<",global_id="<<global_id<<",local_id="<<local_id<<",x="<<x<<",y="<<y<<",z="<<z<<"\n";
00169 #endif
00170         for( int j = 1; j < np_; ++j ) {
00171           (*B_bar_)[j][i] = std::sin(j*pi*z/L);
00172 #ifdef GLPAPP_ADVDIFFREACT_OPTMODEL_DUMP_STUFF
00173           *out << "\nB("<<i<<","<<j<<")="<<(*B_bar_)[j][i]<<"\n";
00174 #endif
00175         }
00176       }
00177       if(normalizeBasis) {
00178 #ifdef HAVE_THYRA_EPETRAEXT
00179         //
00180         // Use modified Gram-Schmidt to create an orthonormal version of B_bar!
00181         //
00182         Teuchos::RefCountPtr<Thyra::MultiVectorBase<double> >
00183           thyra_B_bar = Thyra::create_MultiVector(
00184             B_bar_
00185             ,Thyra::create_VectorSpace(Teuchos::rcp(new Epetra_Map(*map_p_bar_)))
00186             ,Thyra::create_VectorSpace(Teuchos::rcp(new Epetra_Map(*map_p_[p_bndy_idx])))
00187             ),
00188           thyra_fact_R;
00189         Thyra::sillyModifiedGramSchmidt(thyra_B_bar.ptr(), Teuchos::outArg(thyra_fact_R));
00190         Thyra::scale(double(numBndyNodes)/double(np_),  thyra_B_bar.ptr()); // Each row should sum to around one!
00191         // We just discard the "R" factory thyra_fact_R ...
00192 #else // HAVE_THYRA_EPETRAEXT
00193         TEUCHOS_TEST_FOR_EXCEPTION(
00194           true,std::logic_error
00195           ,"Error, can not normalize basis since we do not have Thyra support enabled!"
00196           );
00197 #endif // HAVE_EPETRAEXT_THYRA
00198       }
00199     }
00200 #ifdef GLPAPP_ADVDIFFREACT_OPTMODEL_DUMP_STUFF
00201     *out << "\nB_bar =\n\n";
00202     B_bar_->Print(Teuchos::OSTab(out).o());
00203 #endif
00204   }
00205   else {
00206     // B_bar = I implicitly!
00207     map_p_[p_bndy_idx] = map_p_bar_;
00208   }
00209   //
00210   // Create vectors
00211   //
00212   x0_ = Teuchos::rcp(new Epetra_Vector(*map_x_));
00213   //xL_ = Teuchos::rcp(new Epetra_Vector(*map_x_));
00214   //xU_ = Teuchos::rcp(new Epetra_Vector(*map_x_));
00215   p0_.resize(Np_);
00216   p0_[p_bndy_idx] = Teuchos::rcp(new Epetra_Vector(*map_p_[p_bndy_idx]));
00217   p0_[p_rx_idx] = Teuchos::rcp(new Epetra_Vector(*map_p_[p_rx_idx]));
00218   pL_.resize(Np_);
00219   pU_.resize(Np_);
00220   //pL_ = Teuchos::rcp(new Epetra_Vector(*map_p_));
00221   //pU_ = Teuchos::rcp(new Epetra_Vector(*map_p_));
00222   //
00223   // Initialize the vectors
00224   //
00225   x0_->PutScalar(x0);
00226   p0_[p_bndy_idx]->PutScalar(p0);
00227   p0_[p_rx_idx]->PutScalar(reactionRate);
00228   //
00229   // Initialize the graph for W
00230   //
00231   dat_->computeNpy(x0_);
00232   //if(dumpAll) { *out << "\nA =\n"; { Teuchos::OSTab tab(out); dat_->getA()->Print(*out); } }
00233   //if(dumpAll) { *out << "\nNpy =\n"; {  Teuchos::OSTab tab(out); dat_->getNpy()->Print(*out); } }
00234   W_graph_ = Teuchos::rcp(new Epetra_CrsGraph(dat_->getA()->Graph())); // Assume A and Npy have same graph!
00235   //
00236   // Get default objective matching vector q
00237   //
00238   q_ = Teuchos::rcp(new Epetra_Vector(*(*dat_->getq())(0))); // From Epetra_FEVector to Epetra_Vector!
00239   //
00240   isInitialized_ = true;
00241 #ifdef GLPAPP_ADVDIFFREACT_OPTMODEL_DUMP_STUFF
00242   *out << "\nLeaving AdvDiffReactOptModel::AdvDiffReactOptModel(...) ...\n";
00243 #endif
00244 }
00245 
00246 void AdvDiffReactOptModel::set_q( Teuchos::RefCountPtr<const Epetra_Vector> const& q )
00247 {
00248   q_ = q;
00249 #ifdef GLPAPP_ADVDIFFREACT_OPTMODEL_DUMP_STUFF
00250   Teuchos::RefCountPtr<Teuchos::FancyOStream>
00251     dout = Teuchos::VerboseObjectBase::getDefaultOStream();
00252   Teuchos::OSTab dtab(dout);
00253   *dout << "\nIn AdvDiffReactOptModel::set_q(q):  q =\n\n";
00254   q_->Print(Teuchos::OSTab(dout).o());
00255 #endif
00256 }
00257 
00258 Teuchos::RefCountPtr<GLpApp::GLpYUEpetraDataPool>
00259 AdvDiffReactOptModel::getDataPool()
00260 {
00261   return dat_;
00262 }
00263 
00264 Teuchos::RefCountPtr<const Epetra_MultiVector>
00265 AdvDiffReactOptModel::get_B_bar() const
00266 {
00267   return B_bar_;
00268 }
00269 
00270 // Overridden from EpetraExt::ModelEvaluator
00271 
00272 Teuchos::RefCountPtr<const Epetra_Map>
00273 AdvDiffReactOptModel::get_x_map() const
00274 {
00275   return map_x_;
00276 }
00277 
00278 Teuchos::RefCountPtr<const Epetra_Map>
00279 AdvDiffReactOptModel::get_f_map() const
00280 {
00281   return map_x_;
00282 }
00283 
00284 Teuchos::RefCountPtr<const Epetra_Map>
00285 AdvDiffReactOptModel::get_p_map(int l) const
00286 {
00287   TEUCHOS_TEST_FOR_EXCEPT(!(0<=l<=Np_));
00288   return map_p_[l];
00289 }
00290 
00291 Teuchos::RefCountPtr<const Epetra_Map>
00292 AdvDiffReactOptModel::get_g_map(int j) const
00293 {
00294   TEUCHOS_TEST_FOR_EXCEPT(j!=0);
00295   return map_g_;
00296 }
00297 
00298 Teuchos::RefCountPtr<const Epetra_Vector>
00299 AdvDiffReactOptModel::get_x_init() const
00300 {
00301   return x0_;
00302 }
00303 
00304 Teuchos::RefCountPtr<const Epetra_Vector>
00305 AdvDiffReactOptModel::get_p_init(int l) const
00306 {
00307   TEUCHOS_TEST_FOR_EXCEPT(!(0<=l<=Np_));
00308   return p0_[l];
00309 }
00310 
00311 Teuchos::RefCountPtr<const Epetra_Vector>
00312 AdvDiffReactOptModel::get_x_lower_bounds() const
00313 {
00314   return xL_;
00315 }
00316 
00317 Teuchos::RefCountPtr<const Epetra_Vector>
00318 AdvDiffReactOptModel::get_x_upper_bounds() const
00319 {
00320   return xU_;
00321 }
00322 
00323 Teuchos::RefCountPtr<const Epetra_Vector>
00324 AdvDiffReactOptModel::get_p_lower_bounds(int l) const
00325 {
00326   TEUCHOS_TEST_FOR_EXCEPT(!(0<=l<=Np_));
00327   return pL_[l];
00328 }
00329 
00330 Teuchos::RefCountPtr<const Epetra_Vector>
00331 AdvDiffReactOptModel::get_p_upper_bounds(int l) const
00332 {
00333   TEUCHOS_TEST_FOR_EXCEPT(!(0<=l<=Np_));
00334   return pU_[l];
00335 }
00336 
00337 Teuchos::RefCountPtr<Epetra_Operator>
00338 AdvDiffReactOptModel::create_W() const
00339 {
00340   return Teuchos::rcp(new Epetra_CrsMatrix(::Copy,*W_graph_));
00341 }
00342 
00343 Teuchos::RefCountPtr<Epetra_Operator>
00344 AdvDiffReactOptModel::create_DfDp_op(int l) const
00345 {
00346   TEUCHOS_TEST_FOR_EXCEPT(l!=0);
00347   return Teuchos::rcp(new Epetra_CrsMatrix(::Copy,dat_->getB()->Graph()));
00348   // See DfDp in evalModel(...) below for details
00349 }
00350 
00351 EpetraExt::ModelEvaluator::InArgs
00352 AdvDiffReactOptModel::createInArgs() const
00353 {
00354   InArgsSetup inArgs;
00355   inArgs.setModelEvalDescription(this->description());
00356   inArgs.set_Np(2);
00357   inArgs.setSupports(IN_ARG_x,true);
00358   return inArgs;
00359 }
00360 
00361 EpetraExt::ModelEvaluator::OutArgs
00362 AdvDiffReactOptModel::createOutArgs() const
00363 {
00364   OutArgsSetup outArgs;
00365   outArgs.setModelEvalDescription(this->description());
00366   outArgs.set_Np_Ng(2,1);
00367   outArgs.setSupports(OUT_ARG_f,true);
00368   outArgs.setSupports(OUT_ARG_W,true);
00369   outArgs.set_W_properties(
00370     DerivativeProperties(
00371       DERIV_LINEARITY_NONCONST
00372       ,DERIV_RANK_FULL
00373       ,true // supportsAdjoint
00374       )
00375     );
00376   if(supportDerivatives_) {
00377     outArgs.setSupports(
00378       OUT_ARG_DfDp,0
00379       ,( np_ > 0
00380          ? DerivativeSupport(DERIV_MV_BY_COL)
00381          : DerivativeSupport(DERIV_LINEAR_OP,DERIV_MV_BY_COL)
00382         )
00383       );
00384     outArgs.set_DfDp_properties(
00385       0,DerivativeProperties(
00386         DERIV_LINEARITY_CONST
00387         ,DERIV_RANK_DEFICIENT
00388         ,true // supportsAdjoint
00389         )
00390       );
00391     outArgs.setSupports(OUT_ARG_DgDx,0,DERIV_TRANS_MV_BY_ROW);
00392     outArgs.set_DgDx_properties(
00393       0,DerivativeProperties(
00394         DERIV_LINEARITY_NONCONST
00395         ,DERIV_RANK_DEFICIENT
00396         ,true // supportsAdjoint
00397         )
00398       );
00399     outArgs.setSupports(OUT_ARG_DgDp,0,0,DERIV_TRANS_MV_BY_ROW);
00400     outArgs.set_DgDp_properties(
00401       0,0,DerivativeProperties(
00402         DERIV_LINEARITY_NONCONST
00403         ,DERIV_RANK_DEFICIENT
00404         ,true // supportsAdjoint
00405         )
00406       );
00407   }
00408   return outArgs;
00409 }
00410 
00411 void AdvDiffReactOptModel::evalModel( const InArgs& inArgs, const OutArgs& outArgs ) const
00412 {
00413   Teuchos::TimeMonitor evalTimerMonitor(*evalTimer);
00414 #ifdef GLPAPP_ADVDIFFREACT_OPTMODEL_DUMP_STUFF
00415   Teuchos::RefCountPtr<Teuchos::FancyOStream>
00416     dout = Teuchos::VerboseObjectBase::getDefaultOStream();
00417   Teuchos::OSTab dtab(dout);
00418   *dout << "\nEntering AdvDiffReactOptModel::evalModel(...) ...\n";
00419 #endif
00420   using Teuchos::dyn_cast;
00421   using Teuchos::rcp_dynamic_cast;
00422   //
00423   Teuchos::RefCountPtr<Teuchos::FancyOStream> out = this->getOStream();
00424   const bool trace = ( static_cast<int>(this->getVerbLevel()) >= static_cast<int>(Teuchos::VERB_LOW) );
00425   const bool dumpAll = ( static_cast<int>(this->getVerbLevel()) >= static_cast<int>(Teuchos::VERB_EXTREME) );
00426   //
00427   Teuchos::OSTab tab(out);
00428   if(out.get() && trace) *out << "\n*** Entering AdvDiffReactOptModel::evalModel(...) ...\n"; 
00429   //
00430   // Get the input arguments
00431   //
00432   const Epetra_Vector *p_in = inArgs.get_p(p_bndy_idx).get();
00433   const Epetra_Vector &p = (p_in ? *p_in : *p0_[p_bndy_idx]);
00434   const Epetra_Vector *rx_vec_in = inArgs.get_p(p_rx_idx).get();
00435   const double        reactionRate = ( rx_vec_in ? (*rx_vec_in)[0] : (*p0_[p_rx_idx])[0] );
00436   const Epetra_Vector &x = *inArgs.get_x();
00437 #ifdef GLPAPP_ADVDIFFREACT_OPTMODEL_DUMP_STUFF
00438     *dout << "\nx =\n\n";
00439     x.Print(Teuchos::OSTab(dout).o());
00440     *dout << "\np =\n\n";
00441     p.Print(Teuchos::OSTab(dout).o());
00442 #endif
00443   //
00444   // Get the output arguments
00445   //
00446   Epetra_Vector       *f_out = outArgs.get_f().get();
00447   Epetra_Vector       *g_out = outArgs.get_g(0).get();
00448   Epetra_Operator     *W_out = outArgs.get_W().get();
00449   Derivative          DfDp_out;
00450   Epetra_MultiVector  *DgDx_trans_out = 0;
00451   Epetra_MultiVector  *DgDp_trans_out = 0;
00452   if(supportDerivatives_) {
00453     DfDp_out = outArgs.get_DfDp(0);
00454     DgDx_trans_out = get_DgDx_mv(0,outArgs,DERIV_TRANS_MV_BY_ROW).get();
00455     DgDp_trans_out = get_DgDp_mv(0,0,outArgs,DERIV_TRANS_MV_BY_ROW).get();
00456   }
00457   //
00458   // Precompute some shared quantities
00459   //
00460   // p_bar = B_bar*p
00461   Teuchos::RefCountPtr<const Epetra_Vector> p_bar;
00462   if(np_ > 0) {
00463     Teuchos::TimeMonitor p_bar_TimerMonitor(*p_bar_Timer);
00464     Teuchos::RefCountPtr<Epetra_Vector> _p_bar;
00465     _p_bar = Teuchos::rcp(new Epetra_Vector(*map_p_bar_));
00466     _p_bar->Multiply('N','N',1.0,*B_bar_,p,0.0);
00467     p_bar = _p_bar;
00468   }
00469   else {
00470     p_bar = Teuchos::rcp(&p,false);
00471   }
00472   // R_p_bar = R*p_bar = R*(B_bar*p)
00473   Teuchos::RefCountPtr<const Epetra_Vector> R_p_bar;
00474   if( g_out || DgDp_trans_out ) {
00475     Teuchos::TimeMonitor R_p_bar_TimerMonitor(*R_p_bar_Timer);
00476       Teuchos::RefCountPtr<Epetra_Vector>
00477       _R_p_bar = Teuchos::rcp(new Epetra_Vector(*map_p_bar_));
00478     dat_->getR()->Multiply(false,*p_bar,*_R_p_bar);
00479     R_p_bar = _R_p_bar;
00480   }
00481   //
00482   // Compute the functions
00483   //
00484   if(f_out) {
00485     //
00486     // f = A*x + reationRate*Ny(x) + B*(B_bar*p)
00487     //
00488     Teuchos::TimeMonitor f_TimerMonitor(*f_Timer);
00489     Epetra_Vector &f = *f_out;
00490     Epetra_Vector Ax(*map_f_);
00491     dat_->getA()->Multiply(false,x,Ax);
00492     f.Update(1.0,Ax,0.0);
00493     if(reactionRate!=0.0) {
00494       dat_->computeNy(Teuchos::rcp(&x,false));
00495       f.Update(reactionRate,*dat_->getNy(),1.0);
00496     }
00497     Epetra_Vector Bp(*map_f_);
00498     dat_->getB()->Multiply(false,*p_bar,Bp);
00499     f.Update(1.0,Bp,-1.0,*dat_->getb(),1.0);
00500   }
00501   if(g_out) {
00502     //
00503     // g = 0.5 * (x-q)'*H*(x-q) + 0.5*regBeta*(B_bar*p)'*R*(B_bar*p)
00504     //
00505     Teuchos::TimeMonitor g_TimerMonitor(*g_Timer);
00506     Epetra_Vector &g = *g_out;
00507     Epetra_Vector xq(x);
00508     xq.Update(-1.0, *q_, 1.0);
00509     Epetra_Vector Hxq(x);
00510     dat_->getH()->Multiply(false,xq,Hxq);
00511     g[0] = 0.5*dot(xq,Hxq) + 0.5*dat_->getbeta()*dot(*p_bar,*R_p_bar);
00512 #ifdef GLPAPP_ADVDIFFREACT_OPTMODEL_DUMP_STUFF
00513     *dout << "q =\n\n";
00514     q_->Print(Teuchos::OSTab(dout).o());
00515     *dout << "x =\n\n";
00516     x.Print(Teuchos::OSTab(dout).o());
00517     *dout << "x-q =\n\n";
00518     xq.Print(Teuchos::OSTab(dout).o());
00519     *dout << "H*(x-q) =\n\n";
00520     Hxq.Print(Teuchos::OSTab(dout).o());
00521     *dout << "0.5*dot(xq,Hxq) = " << (0.5*dot(xq,Hxq)) << "\n";
00522     *dout << "0.5*regBeta*(B_bar*p)'*R*(B_bar*p) = " << (0.5*dat_->getbeta()*dot(*p_bar,*R_p_bar)) << "\n";
00523 #endif
00524   }
00525   if(W_out) {
00526     //
00527     // W = A + reationRate*Npy(x)
00528     //
00529     Teuchos::TimeMonitor W_TimerMonitor(*W_Timer);
00530     Epetra_CrsMatrix &DfDx = dyn_cast<Epetra_CrsMatrix>(*W_out);
00531     if(reactionRate!=0.0)
00532       dat_->computeNpy(Teuchos::rcp(&x,false));
00533     Teuchos::RefCountPtr<Epetra_CrsMatrix>
00534       dat_A = dat_->getA(),
00535       dat_Npy = dat_->getNpy();
00536     const int numMyRows = dat_A->NumMyRows();
00537     for( int i = 0; i < numMyRows; ++i ) {
00538       int dat_A_num_row_entries=0; double *dat_A_row_vals=0; int *dat_A_row_inds=0;
00539       dat_A->ExtractMyRowView(i,dat_A_num_row_entries,dat_A_row_vals,dat_A_row_inds);
00540       int DfDx_num_row_entries=0; double *DfDx_row_vals=0; int *DfDx_row_inds=0;
00541       DfDx.ExtractMyRowView(i,DfDx_num_row_entries,DfDx_row_vals,DfDx_row_inds);
00542 #ifdef TEUCHOS_DEBUG
00543       TEUCHOS_TEST_FOR_EXCEPT(DfDx_num_row_entries!=dat_A_num_row_entries);
00544 #endif
00545       if(reactionRate!=0.0) {
00546         int dat_Npy_num_row_entries=0; double *dat_Npy_row_vals=0; int *dat_Npy_row_inds=0;
00547         dat_Npy->ExtractMyRowView(i,dat_Npy_num_row_entries,dat_Npy_row_vals,dat_Npy_row_inds);
00548 #ifdef TEUCHOS_DEBUG
00549         TEUCHOS_TEST_FOR_EXCEPT(dat_A_num_row_entries!=dat_Npy_num_row_entries);
00550 #endif
00551         for(int k = 0; k < DfDx_num_row_entries; ++k) {
00552 #ifdef TEUCHOS_DEBUG
00553           TEUCHOS_TEST_FOR_EXCEPT(dat_A_row_inds[k]!=dat_Npy_row_inds[k]||dat_A_row_inds[k]!=DfDx_row_inds[k]);
00554 #endif
00555           DfDx_row_vals[k] = dat_A_row_vals[k] + reactionRate * dat_Npy_row_vals[k];
00556         }
00557       }
00558       else {
00559         for(int k = 0; k < DfDx_num_row_entries; ++k) {
00560 #ifdef TEUCHOS_DEBUG
00561           TEUCHOS_TEST_FOR_EXCEPT(dat_A_row_inds[k]!=DfDx_row_inds[k]);
00562 #endif
00563           DfDx_row_vals[k] = dat_A_row_vals[k];
00564         }
00565       }
00566     }
00567   }
00568   if(!DfDp_out.isEmpty()) {
00569     if(out.get() && trace) *out << "\nComputing DfDp ...\n"; 
00570     //
00571     // DfDp = B*B_bar
00572     //
00573     Teuchos::TimeMonitor DfDp_TimerMonitor(*DfDp_Timer);
00574     Epetra_CrsMatrix   *DfDp_op = NULL;
00575     Epetra_MultiVector *DfDp_mv = NULL;
00576     if(out.get() && dumpAll)
00577     { *out << "\nB =\n"; { Teuchos::OSTab tab(out); dat_->getB()->Print(*out); } }
00578     if(DfDp_out.getLinearOp().get()) {
00579       DfDp_op = &dyn_cast<Epetra_CrsMatrix>(*DfDp_out.getLinearOp());
00580     }
00581     else {
00582       DfDp_mv = &*DfDp_out.getDerivativeMultiVector().getMultiVector();
00583       DfDp_mv->PutScalar(0.0);
00584     }
00585     Teuchos::RefCountPtr<Epetra_CrsMatrix>
00586       dat_B = dat_->getB();
00587     if(np_ > 0) {
00588       //
00589       // We only support a Multi-vector form when we have a non-idenity basis
00590       // matrix B_bar for p!
00591       //
00592       TEUCHOS_TEST_FOR_EXCEPT(DfDp_mv==NULL);
00593       dat_B->Multiply(false,*B_bar_,*DfDp_mv);
00594     }
00595     else {
00596       //
00597       // Note: We copy from B every time in order to be safe.  Note that since
00598       // the client knows that B is constant (sense we told them so in
00599       // createOutArgs()) then it should only compute this matrix once and keep
00600       // it if it is smart.
00601       //
00602       // Note: We support both the CrsMatrix and MultiVector form of this object
00603       // to make things easier for the client.
00604       //
00605       if(DfDp_op) {
00606         const int numMyRows = dat_B->NumMyRows();
00607         for( int i = 0; i < numMyRows; ++i ) {
00608           int dat_B_num_row_entries=0; double *dat_B_row_vals=0; int *dat_B_row_inds=0;
00609           dat_B->ExtractMyRowView(i,dat_B_num_row_entries,dat_B_row_vals,dat_B_row_inds);
00610           int DfDp_num_row_entries=0; double *DfDp_row_vals=0; int *DfDp_row_inds=0;
00611           DfDp_op->ExtractMyRowView(i,DfDp_num_row_entries,DfDp_row_vals,DfDp_row_inds);
00612 #ifdef TEUCHOS_DEBUG
00613           TEUCHOS_TEST_FOR_EXCEPT(DfDp_num_row_entries!=dat_B_num_row_entries);
00614 #endif
00615           for(int k = 0; k < DfDp_num_row_entries; ++k) {
00616 #ifdef TEUCHOS_DEBUG
00617             TEUCHOS_TEST_FOR_EXCEPT(dat_B_row_inds[k]!=DfDp_row_inds[k]);
00618 #endif
00619             DfDp_row_vals[k] = dat_B_row_vals[k];
00620           }
00621           // ToDo: The above code should be put in a utility function called copyValues(...)!
00622         }
00623       }
00624       else if(DfDp_mv) {
00625         // We must do a mat-vec to get this since we can't just copy out the
00626         // matrix entries since the domain map may be different from the
00627         // column map!  I learned this the very very hard way!  I am using
00628         // Thyra wrappers here since I can't figure out for the life of me how
00629         // to do this cleanly with Epetra alone!
00630         Teuchos::RefCountPtr<Epetra_Vector>
00631           etaVec = Teuchos::rcp(new Epetra_Vector(*map_p_bar_));
00632         Teuchos::RefCountPtr<const Thyra::VectorSpaceBase<double> >
00633           space_p_bar = Thyra::create_VectorSpace(Teuchos::rcp(new Epetra_Map(*map_p_bar_)));
00634         Teuchos::RefCountPtr<Thyra::VectorBase<double> >
00635           thyra_etaVec = Thyra::create_Vector(etaVec,space_p_bar);
00636         for( int i = 0; i < map_p_bar_->NumGlobalElements(); ++i ) {
00637           Thyra::assign(thyra_etaVec.ptr(), 0.0);
00638           Thyra::set_ele(i, 1.0, thyra_etaVec.ptr());
00639           dat_B->Multiply(false, *etaVec, *(*DfDp_mv)(i));
00640         };
00641       }
00642     }
00643 #ifdef GLPAPP_ADVDIFFREACT_OPTMODEL_DUMP_STUFF
00644     if(DfDp_op) {
00645       *dout << "\nDfDp_op =\n\n";
00646       DfDp_op->Print(Teuchos::OSTab(dout).o());
00647     }
00648     if(DfDp_mv) {
00649       *dout << "\nDfDp_mv =\n\n";
00650       DfDp_mv->Print(Teuchos::OSTab(dout).o());
00651     }
00652 #endif
00653   }
00654   if(DgDx_trans_out) {
00655     //
00656     // DgDx' = H*(x-q)
00657     //
00658     Teuchos::TimeMonitor DgDx_TimerMonitor(*DgDx_Timer);
00659     Epetra_Vector &DgDx_trans = *(*DgDx_trans_out)(0);
00660     Epetra_Vector xq(x);
00661     xq.Update(-1.0,*q_,1.0);
00662     dat_->getH()->Multiply(false,xq,DgDx_trans);
00663 #ifdef GLPAPP_ADVDIFFREACT_OPTMODEL_DUMP_STUFF
00664     *dout << "q =\n\n";
00665     q_->Print(Teuchos::OSTab(dout).o());
00666     *dout << "x =\n\n";
00667     x.Print(Teuchos::OSTab(dout).o());
00668     *dout << "x-q =\n\n";
00669     xq.Print(Teuchos::OSTab(dout).o());
00670     *dout << "DgDx_trans = H*(x-q) =\n\n";
00671     DgDx_trans.Print(Teuchos::OSTab(dout).o());
00672 #endif
00673   }
00674   if(DgDp_trans_out) {
00675     //
00676     // DgDp' = regBeta*B_bar'*(R*(B_bar*p))
00677     //
00678     Teuchos::TimeMonitor DgDp_TimerMonitor(*DgDp_Timer);
00679     Epetra_Vector &DgDp_trans = *(*DgDp_trans_out)(0);
00680     if(np_ > 0) {
00681       DgDp_trans.Multiply('T','N',dat_->getbeta(),*B_bar_,*R_p_bar,0.0);
00682     }
00683     else {
00684       DgDp_trans.Update(dat_->getbeta(),*R_p_bar,0.0);
00685     }
00686 #ifdef GLPAPP_ADVDIFFREACT_OPTMODEL_DUMP_STUFF
00687     *dout << "\nR_p_bar =\n\n";
00688     R_p_bar->Print(Teuchos::OSTab(dout).o());
00689     if(B_bar_.get()) {
00690       *dout << "\nB_bar =\n\n";
00691       B_bar_->Print(Teuchos::OSTab(dout).o());
00692     }
00693     *dout << "\nDgDp_trans =\n\n";
00694     DgDp_trans.Print(Teuchos::OSTab(dout).o());
00695 #endif
00696   }
00697   if(out.get() && trace) *out << "\n*** Leaving AdvDiffReactOptModel::evalModel(...) ...\n"; 
00698 #ifdef GLPAPP_ADVDIFFREACT_OPTMODEL_DUMP_STUFF
00699   *dout << "\nLeaving AdvDiffReactOptModel::evalModel(...) ...\n";
00700 #endif
00701 }
00702 
00703 } // namespace GLpApp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines