GLpApp_AdvDiffReactOptModel.cpp

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

Generated on Wed May 12 21:40:37 2010 for EpetraExt by  doxygen 1.4.7