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
00009
00010
00011
00012
00013 #include "Thyra_EpetraThyraWrappers.hpp"
00014 #include "Thyra_VectorStdOps.hpp"
00015
00016 #ifdef HAVE_THYRA_EPETRAEXT
00017
00018 #include "sillyModifiedGramSchmidt.hpp"
00019 #include "Thyra_MultiVectorStdOps.hpp"
00020 #include "Thyra_SpmdVectorSpaceBase.hpp"
00021 #endif // HAVE_THYRA_EPETRAEXT
00022
00023
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 }
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
00056 ,const double len_y
00057 ,const int local_nx
00058 ,const int local_ny
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
00078
00079 dat_ = Teuchos::rcp(new GLpApp::GLpYUEpetraDataPool(comm,beta,len_x,len_y,local_nx,local_ny,meshFile,false));
00080
00081
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
00091
00092 if( np_ > 0 ) {
00093
00094
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);
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
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);
00146
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
00162 map_p_[p_bndy_idx] = map_p_bar_;
00163 }
00164
00165
00166
00167 x0_ = Teuchos::rcp(new Epetra_Vector(*map_x_));
00168
00169
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
00176
00177
00178
00179
00180 x0_->PutScalar(x0);
00181 p0_[p_bndy_idx]->PutScalar(p0);
00182 p0_[p_rx_idx]->PutScalar(reactionRate);
00183
00184
00185
00186 dat_->computeNpy(x0_);
00187
00188
00189 W_graph_ = Teuchos::rcp(new Epetra_CrsGraph(dat_->getA()->Graph()));
00190
00191
00192
00193 q_ = Teuchos::rcp(new Epetra_Vector(*(*dat_->getq())(0)));
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
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
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
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
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
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
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
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
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
00414
00415
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
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
00438
00439 if(f_out) {
00440
00441
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
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
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
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
00545
00546
00547 TEST_FOR_EXCEPT(DfDp_mv==NULL);
00548 dat_B->Multiply(false,*B_bar_,*DfDp_mv);
00549 }
00550 else {
00551
00552
00553
00554
00555
00556
00557
00558
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
00577 }
00578 }
00579 else if(DfDp_mv) {
00580
00581
00582
00583
00584
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
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
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 }