47 #include "Teuchos_ScalarTraits.hpp" 48 #include "Teuchos_VerboseObject.hpp" 49 #include "Teuchos_TimeMonitor.hpp" 56 #include "Thyra_EpetraThyraWrappers.hpp" 57 #include "Thyra_VectorStdOps.hpp" 59 #ifdef HAVE_THYRA_EPETRAEXT 61 #include "sillyModifiedGramSchmidt.hpp" 62 #include "Thyra_MultiVectorStdOps.hpp" 63 #include "Thyra_SpmdVectorSpaceBase.hpp" 64 #endif // HAVE_THYRA_EPETRAEXT 70 Teuchos::RCP<Teuchos::Time>
71 initalizationTimer = Teuchos::TimeMonitor::getNewTimer(
"AdvDiffReact:Init"),
72 evalTimer = Teuchos::TimeMonitor::getNewTimer(
"AdvDiffReact:Eval"),
73 p_bar_Timer = Teuchos::TimeMonitor::getNewTimer(
"AdvDiffReact:Eval:p_bar"),
74 R_p_bar_Timer = Teuchos::TimeMonitor::getNewTimer(
"AdvDiffReact:Eval:R_p_bar"),
75 f_Timer = Teuchos::TimeMonitor::getNewTimer(
"AdvDiffReact:Eval:f"),
76 g_Timer = Teuchos::TimeMonitor::getNewTimer(
"AdvDiffReact:Eval:g"),
77 W_Timer = Teuchos::TimeMonitor::getNewTimer(
"AdvDiffReact:Eval:W"),
78 DfDp_Timer = Teuchos::TimeMonitor::getNewTimer(
"AdvDiffReact:Eval:DfDp"),
79 DgDx_Timer = Teuchos::TimeMonitor::getNewTimer(
"AdvDiffReact:Eval:DgDx"),
80 DgDp_Timer = Teuchos::TimeMonitor::getNewTimer(
"AdvDiffReact:Eval:DgDp");
82 inline double sqr(
const double& s ) {
return s*s; }
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());
192 #else // HAVE_THYRA_EPETRAEXT 193 TEUCHOS_TEST_FOR_EXCEPTION(
194 true,std::logic_error
195 ,
"Error, can not normalize basis since we do not have Thyra support enabled!" 197 #endif // HAVE_EPETRAEXT_THYRA 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";
249 #ifdef GLPAPP_ADVDIFFREACT_OPTMODEL_DUMP_STUFF 250 Teuchos::RCP<Teuchos::FancyOStream>
251 dout = Teuchos::VerboseObjectBase::getDefaultOStream();
252 Teuchos::OSTab dtab(dout);
253 *dout <<
"\nIn AdvDiffReactOptModel::set_q(q): q =\n\n";
254 q_->Print(Teuchos::OSTab(dout).o());
258 Teuchos::RCP<GLpApp::GLpYUEpetraDataPool>
264 Teuchos::RCP<const Epetra_MultiVector>
272 Teuchos::RCP<const Epetra_Map>
278 Teuchos::RCP<const Epetra_Map>
284 Teuchos::RCP<const Epetra_Map>
287 TEUCHOS_TEST_FOR_EXCEPT(!(0<=l<=
Np_));
291 Teuchos::RCP<const Epetra_Map>
294 TEUCHOS_TEST_FOR_EXCEPT(j!=0);
298 Teuchos::RCP<const Epetra_Vector>
304 Teuchos::RCP<const Epetra_Vector>
307 TEUCHOS_TEST_FOR_EXCEPT(!(0<=l<=
Np_));
311 Teuchos::RCP<const Epetra_Vector>
317 Teuchos::RCP<const Epetra_Vector>
323 Teuchos::RCP<const Epetra_Vector>
326 TEUCHOS_TEST_FOR_EXCEPT(!(0<=l<=
Np_));
330 Teuchos::RCP<const Epetra_Vector>
333 TEUCHOS_TEST_FOR_EXCEPT(!(0<=l<=
Np_));
337 Teuchos::RCP<Epetra_Operator>
343 Teuchos::RCP<Epetra_Operator>
346 TEUCHOS_TEST_FOR_EXCEPT(l!=0);
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());
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);
493 if(reactionRate!=0.0) {
494 dat_->computeNy(Teuchos::rcp(&x,
false));
498 dat_->getB()->Multiply(
false,*p_bar,Bp);
505 Teuchos::TimeMonitor g_TimerMonitor(*g_Timer);
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;
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];
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); } }
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";
Teuchos::RCP< Epetra_Vector > x0_
RCP_Eptra_Vector_Array_t pU_
int Dot(const Epetra_MultiVector &A, double *Result) const
Teuchos::RCP< GLpApp::GLpYUEpetraDataPool > getDataPool()
void set_W_properties(const DerivativeProperties &properties)
Teuchos::RCP< const Epetra_Vector > get_x_lower_bounds() const
AdvDiffReactOptModel(const Teuchos::RCP< const Epetra_Comm > &comm, const double beta, const double len_x, const double len_y, const int local_nx, const int local_ny, const char meshFile[], const int np, const double x0, const double p0, const double reactionRate, const bool normalizeBasis, const bool supportDerivatives)
Constructor.
Teuchos::RCP< Epetra_CrsGraph > W_graph_
Teuchos::RCP< const Epetra_Map > map_f_
void setSupports(EOutArgsMembers arg, bool supports=true)
Teuchos::RCP< Epetra_Operator > create_DfDp_op(int l) const
Teuchos::RCP< const Epetra_Map > get_x_map() const
virtual void Print(std::ostream &os) const
Teuchos::RCP< Epetra_MultiVector > B_bar_
Teuchos::RCP< const Epetra_Vector > get_x_upper_bounds() const
Teuchos::RCP< Epetra_MultiVector > get_DgDp_mv(const int j, const int l, const ModelEvaluator::OutArgs &outArgs, const ModelEvaluator::EDerivativeMultiVectorOrientation mvOrientation)
int PutScalar(double ScalarConstant)
Teuchos::RCP< const Epetra_Vector > get_p_lower_bounds(int l) const
void setSupports(EInArgsMembers arg, bool supports=true)
Teuchos::RCP< GLpApp::GLpYUEpetraDataPool > dat_
InArgs createInArgs() const
Simple aggregate class that stores a derivative object as a general linear operator or as a multi-vec...
static const int p_rx_idx
int ExtractMyRowView(int MyRow, int &NumEntries, double *&Values, int *&Indices) const
Teuchos::RCP< const Epetra_Map > map_g_
Teuchos::RCP< Epetra_Operator > create_W() const
Teuchos::RCP< const Epetra_Vector > get_x() const
Set solution vector Taylor polynomial.
Teuchos::RCP< const Epetra_Map > get_g_map(int j) const
.
Teuchos::RCP< const Epetra_Map > get_p_map(int l) const
.
void set_DgDp_properties(int j, int l, const DerivativeProperties &properties)
Derivative get_DfDp(int l) const
RCP_Eptra_Map_Array_t map_p_
RCP_Eptra_Vector_Array_t p0_
void set_DfDp_properties(int l, const DerivativeProperties &properties)
void setModelEvalDescription(const std::string &modelEvalDescription)
Teuchos::RCP< Epetra_Operator > get_W() const
Teuchos::RCP< const Epetra_Map > map_x_
Teuchos::RCP< Epetra_Vector > xU_
Teuchos::RCP< Epetra_MultiVector > get_DgDx_mv(const int j, const ModelEvaluator::OutArgs &outArgs, const ModelEvaluator::EDerivativeMultiVectorOrientation mvOrientation)
Teuchos::RCP< Epetra_Vector > xL_
Teuchos::RCP< const Epetra_MultiVector > get_B_bar() const
Teuchos::RCP< const Epetra_Vector > get_p_init(int l) const
void setModelEvalDescription(const std::string &modelEvalDescription)
DerivativeMultiVector getDerivativeMultiVector() const
Teuchos::RCP< const Epetra_Vector > get_p(int l) const
virtual void Print(std::ostream &os) const
Teuchos::RCP< const Epetra_Vector > q_
int Update(double ScalarA, const Epetra_MultiVector &A, double ScalarThis)
int Multiply(char TransA, char TransB, double ScalarAB, const Epetra_MultiVector &A, const Epetra_MultiVector &B, double ScalarThis)
Teuchos::RCP< Epetra_MultiVector > getMultiVector() const
void evalModel(const InArgs &inArgs, const OutArgs &outArgs) const
void set_q(Teuchos::RCP< const Epetra_Vector > const &q)
Evaluation< Epetra_Vector > get_g(int j) const
Get g(j) where 0 <= j && j < this->Ng().
void set_DgDx_properties(int j, const DerivativeProperties &properties)
Teuchos::RCP< const Epetra_Map > get_f_map() const
Evaluation< Epetra_Vector > get_f() const
Teuchos::RCP< Epetra_Operator > getLinearOp() const
Teuchos::RCP< const Epetra_Vector > get_p_upper_bounds(int l) const
OutArgs createOutArgs() const
RCP_Eptra_Vector_Array_t pL_
static const int p_bndy_idx
Teuchos::RCP< const Epetra_Vector > get_x_init() const
void set_Np_Ng(int Np, int Ng)
Teuchos::RCP< const Epetra_Map > map_p_bar_