96 const Teuchos::RCP<const Epetra_Comm> &comm
102 ,
const char meshFile[]
106 ,
const double reactionRate
107 ,
const bool normalizeBasis
108 ,
const bool supportDerivatives
110 : supportDerivatives_(supportDerivatives), np_(np)
112 Teuchos::TimeMonitor initalizationTimerMonitor(*initalizationTimer);
113#ifdef GLPAPP_ADVDIFFREACT_OPTMODEL_DUMP_STUFF
114 Teuchos::RCP<Teuchos::FancyOStream>
115 out = Teuchos::VerboseObjectBase::getDefaultOStream();
116 Teuchos::OSTab tab(out);
117 *out <<
"\nEntering AdvDiffReactOptModel::AdvDiffReactOptModel(...) ...\n";
141 Epetra_Map overlapmap(-1,pindx.M(),
const_cast<int*
>(pindx.A()),1,*comm);
142 const double pi = 2.0 * std::asin(1.0);
143#ifdef GLPAPP_ADVDIFFREACT_OPTMODEL_DUMP_STUFF
144 *out <<
"\npi="<<pi<<
"\n";
148 (*B_bar_)(0)->PutScalar(1.0);
150 const int numBndyNodes =
map_p_bar_->NumMyElements();
151 const int *bndyNodeGlobalIDs =
map_p_bar_->MyGlobalElements();
152 for(
int i = 0; i < numBndyNodes; ++i ) {
153#ifdef GLPAPP_ADVDIFFREACT_OPTMODEL_DUMP_STUFF
154 const int global_id = bndyNodeGlobalIDs[i];
156 const int local_id = overlapmap.LID(bndyNodeGlobalIDs[i]);
157 const double x = ipcoords(local_id,0), y = ipcoords(local_id,1);
158 double z = -1.0, L = -1.0;
159 if( x < 1e-10 || len_x - 1e-10 < x ) {
167#ifdef GLPAPP_ADVDIFFREACT_OPTMODEL_DUMP_STUFF
168 *out <<
"\ni="<<i<<
",global_id="<<global_id<<
",local_id="<<local_id<<
",x="<<x<<
",y="<<y<<
",z="<<z<<
"\n";
170 for(
int j = 1; j <
np_; ++j ) {
171 (*B_bar_)[j][i] = std::sin(j*pi*z/L);
172#ifdef GLPAPP_ADVDIFFREACT_OPTMODEL_DUMP_STUFF
173 *out <<
"\nB("<<i<<
","<<j<<
")="<<(*B_bar_)[j][i]<<
"\n";
178#ifdef HAVE_THYRA_EPETRAEXT
182 Teuchos::RCP<Thyra::MultiVectorBase<double> >
183 thyra_B_bar = Thyra::create_MultiVector(
189 Thyra::sillyModifiedGramSchmidt(thyra_B_bar.ptr(), Teuchos::outArg(thyra_fact_R));
190 Thyra::scale(
double(numBndyNodes)/
double(
np_), thyra_B_bar.ptr());
193 TEUCHOS_TEST_FOR_EXCEPTION(
194 true,std::logic_error
195 ,
"Error, can not normalize basis since we do not have Thyra support enabled!"
200#ifdef GLPAPP_ADVDIFFREACT_OPTMODEL_DUMP_STUFF
201 *out <<
"\nB_bar =\n\n";
202 B_bar_->Print(Teuchos::OSTab(out).o());
241#ifdef GLPAPP_ADVDIFFREACT_OPTMODEL_DUMP_STUFF
242 *out <<
"\nLeaving AdvDiffReactOptModel::AdvDiffReactOptModel(...) ...\n";
413 Teuchos::TimeMonitor evalTimerMonitor(*evalTimer);
414#ifdef GLPAPP_ADVDIFFREACT_OPTMODEL_DUMP_STUFF
415 Teuchos::RCP<Teuchos::FancyOStream>
416 dout = Teuchos::VerboseObjectBase::getDefaultOStream();
417 Teuchos::OSTab dtab(dout);
418 *dout <<
"\nEntering AdvDiffReactOptModel::evalModel(...) ...\n";
420 using Teuchos::dyn_cast;
421 using Teuchos::rcp_dynamic_cast;
423 Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
424 const bool trace = (
static_cast<int>(this->getVerbLevel()) >=
static_cast<int>(Teuchos::VERB_LOW) );
425 const bool dumpAll = (
static_cast<int>(this->getVerbLevel()) >=
static_cast<int>(Teuchos::VERB_EXTREME) );
427 Teuchos::OSTab tab(out);
428 if(out.get() && trace) *out <<
"\n*** Entering AdvDiffReactOptModel::evalModel(...) ...\n";
435 const double reactionRate = ( rx_vec_in ? (*rx_vec_in)[0] : (*
p0_[
p_rx_idx])[0] );
437#ifdef GLPAPP_ADVDIFFREACT_OPTMODEL_DUMP_STUFF
438 *dout <<
"\nx =\n\n";
439 x.
Print(Teuchos::OSTab(dout).o());
440 *dout <<
"\np =\n\n";
441 p.Print(Teuchos::OSTab(dout).o());
453 DfDp_out = outArgs.get_DfDp(0);
461 Teuchos::RCP<const Epetra_Vector> p_bar;
463 Teuchos::TimeMonitor p_bar_TimerMonitor(*p_bar_Timer);
464 Teuchos::RCP<Epetra_Vector> _p_bar;
466 _p_bar->Multiply(
'N',
'N',1.0,*
B_bar_,p,0.0);
470 p_bar = Teuchos::rcp(&p,
false);
473 Teuchos::RCP<const Epetra_Vector> R_p_bar;
474 if( g_out || DgDp_trans_out ) {
475 Teuchos::TimeMonitor R_p_bar_TimerMonitor(*R_p_bar_Timer);
476 Teuchos::RCP<Epetra_Vector>
478 dat_->getR()->Multiply(
false,*p_bar,*_R_p_bar);
488 Teuchos::TimeMonitor f_TimerMonitor(*f_Timer);
491 dat_->getA()->Multiply(
false,x,Ax);
492 f.Update(1.0,Ax,0.0);
493 if(reactionRate!=0.0) {
494 dat_->computeNy(Teuchos::rcp(&x,
false));
495 f.Update(reactionRate,*
dat_->getNy(),1.0);
498 dat_->getB()->Multiply(
false,*p_bar,Bp);
499 f.Update(1.0,Bp,-1.0,*
dat_->getb(),1.0);
505 Teuchos::TimeMonitor g_TimerMonitor(*g_Timer);
508 xq.Update(-1.0, *
q_, 1.0);
510 dat_->getH()->Multiply(
false,xq,Hxq);
511 g[0] = 0.5*dot(xq,Hxq) + 0.5*
dat_->getbeta()*dot(*p_bar,*R_p_bar);
512#ifdef GLPAPP_ADVDIFFREACT_OPTMODEL_DUMP_STUFF
514 q_->Print(Teuchos::OSTab(dout).o());
516 x.Print(Teuchos::OSTab(dout).o());
517 *dout <<
"x-q =\n\n";
518 xq.Print(Teuchos::OSTab(dout).o());
519 *dout <<
"H*(x-q) =\n\n";
520 Hxq.Print(Teuchos::OSTab(dout).o());
521 *dout <<
"0.5*dot(xq,Hxq) = " << (0.5*dot(xq,Hxq)) <<
"\n";
522 *dout <<
"0.5*regBeta*(B_bar*p)'*R*(B_bar*p) = " << (0.5*
dat_->getbeta()*dot(*p_bar,*R_p_bar)) <<
"\n";
529 Teuchos::TimeMonitor W_TimerMonitor(*W_Timer);
531 if(reactionRate!=0.0)
532 dat_->computeNpy(Teuchos::rcp(&x,
false));
533 Teuchos::RCP<Epetra_CrsMatrix>
534 dat_A =
dat_->getA(),
535 dat_Npy =
dat_->getNpy();
536 const int numMyRows = dat_A->NumMyRows();
537 for(
int i = 0; i < numMyRows; ++i ) {
538 int dat_A_num_row_entries=0;
double *dat_A_row_vals=0;
int *dat_A_row_inds=0;
539 dat_A->ExtractMyRowView(i,dat_A_num_row_entries,dat_A_row_vals,dat_A_row_inds);
540 int DfDx_num_row_entries=0;
double *DfDx_row_vals=0;
int *DfDx_row_inds=0;
541 DfDx.ExtractMyRowView(i,DfDx_num_row_entries,DfDx_row_vals,DfDx_row_inds);
543 TEUCHOS_TEST_FOR_EXCEPT(DfDx_num_row_entries!=dat_A_num_row_entries);
545 if(reactionRate!=0.0) {
546 int dat_Npy_num_row_entries=0;
double *dat_Npy_row_vals=0;
int *dat_Npy_row_inds=0;
547 dat_Npy->ExtractMyRowView(i,dat_Npy_num_row_entries,dat_Npy_row_vals,dat_Npy_row_inds);
549 TEUCHOS_TEST_FOR_EXCEPT(dat_A_num_row_entries!=dat_Npy_num_row_entries);
551 for(
int k = 0; k < DfDx_num_row_entries; ++k) {
553 TEUCHOS_TEST_FOR_EXCEPT(dat_A_row_inds[k]!=dat_Npy_row_inds[k]||dat_A_row_inds[k]!=DfDx_row_inds[k]);
555 DfDx_row_vals[k] = dat_A_row_vals[k] + reactionRate * dat_Npy_row_vals[k];
559 for(
int k = 0; k < DfDx_num_row_entries; ++k) {
561 TEUCHOS_TEST_FOR_EXCEPT(dat_A_row_inds[k]!=DfDx_row_inds[k]);
563 DfDx_row_vals[k] = dat_A_row_vals[k];
568 if(!DfDp_out.isEmpty()) {
569 if(out.get() && trace) *out <<
"\nComputing DfDp ...\n";
573 Teuchos::TimeMonitor DfDp_TimerMonitor(*DfDp_Timer);
576 if(out.get() && dumpAll)
577 { *out <<
"\nB =\n"; { Teuchos::OSTab tab(out);
dat_->getB()->Print(*out); } }
578 if(DfDp_out.getLinearOp().get()) {
579 DfDp_op = &dyn_cast<Epetra_CrsMatrix>(*DfDp_out.getLinearOp());
582 DfDp_mv = &*DfDp_out.getDerivativeMultiVector().getMultiVector();
583 DfDp_mv->PutScalar(0.0);
585 Teuchos::RCP<Epetra_CrsMatrix>
586 dat_B =
dat_->getB();
592 TEUCHOS_TEST_FOR_EXCEPT(DfDp_mv==NULL);
593 dat_B->Multiply(
false,*
B_bar_,*DfDp_mv);
606 const int numMyRows = dat_B->NumMyRows();
607 for(
int i = 0; i < numMyRows; ++i ) {
608 int dat_B_num_row_entries=0;
double *dat_B_row_vals=0;
int *dat_B_row_inds=0;
609 dat_B->ExtractMyRowView(i,dat_B_num_row_entries,dat_B_row_vals,dat_B_row_inds);
610 int DfDp_num_row_entries=0;
double *DfDp_row_vals=0;
int *DfDp_row_inds=0;
611 DfDp_op->ExtractMyRowView(i,DfDp_num_row_entries,DfDp_row_vals,DfDp_row_inds);
613 TEUCHOS_TEST_FOR_EXCEPT(DfDp_num_row_entries!=dat_B_num_row_entries);
615 for(
int k = 0; k < DfDp_num_row_entries; ++k) {
617 TEUCHOS_TEST_FOR_EXCEPT(dat_B_row_inds[k]!=DfDp_row_inds[k]);
619 DfDp_row_vals[k] = dat_B_row_vals[k];
630 Teuchos::RCP<Epetra_Vector>
632 Teuchos::RCP<const Thyra::VectorSpaceBase<double> >
634 Teuchos::RCP<Thyra::VectorBase<double> >
635 thyra_etaVec = Thyra::create_Vector(etaVec,space_p_bar);
636 for(
int i = 0; i <
map_p_bar_->NumGlobalElements(); ++i ) {
637 Thyra::assign(thyra_etaVec.ptr(), 0.0);
638 Thyra::set_ele(i, 1.0, thyra_etaVec.ptr());
639 dat_B->Multiply(
false, *etaVec, *(*DfDp_mv)(i));
643#ifdef GLPAPP_ADVDIFFREACT_OPTMODEL_DUMP_STUFF
645 *dout <<
"\nDfDp_op =\n\n";
646 DfDp_op->Print(Teuchos::OSTab(dout).o());
649 *dout <<
"\nDfDp_mv =\n\n";
650 DfDp_mv->Print(Teuchos::OSTab(dout).o());
658 Teuchos::TimeMonitor DgDx_TimerMonitor(*DgDx_Timer);
661 xq.Update(-1.0,*
q_,1.0);
662 dat_->getH()->Multiply(
false,xq,DgDx_trans);
663#ifdef GLPAPP_ADVDIFFREACT_OPTMODEL_DUMP_STUFF
665 q_->Print(Teuchos::OSTab(dout).o());
667 x.Print(Teuchos::OSTab(dout).o());
668 *dout <<
"x-q =\n\n";
669 xq.Print(Teuchos::OSTab(dout).o());
670 *dout <<
"DgDx_trans = H*(x-q) =\n\n";
671 DgDx_trans.Print(Teuchos::OSTab(dout).o());
678 Teuchos::TimeMonitor DgDp_TimerMonitor(*DgDp_Timer);
684 DgDp_trans.Update(
dat_->getbeta(),*R_p_bar,0.0);
686#ifdef GLPAPP_ADVDIFFREACT_OPTMODEL_DUMP_STUFF
687 *dout <<
"\nR_p_bar =\n\n";
688 R_p_bar->Print(Teuchos::OSTab(dout).o());
690 *dout <<
"\nB_bar =\n\n";
691 B_bar_->Print(Teuchos::OSTab(dout).o());
693 *dout <<
"\nDgDp_trans =\n\n";
694 DgDp_trans.Print(Teuchos::OSTab(dout).o());
697 if(out.get() && trace) *out <<
"\n*** Leaving AdvDiffReactOptModel::evalModel(...) ...\n";
698#ifdef GLPAPP_ADVDIFFREACT_OPTMODEL_DUMP_STUFF
699 *dout <<
"\nLeaving AdvDiffReactOptModel::evalModel(...) ...\n";