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"
70Teuchos::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");
82inline 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());
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";
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());
258Teuchos::RCP<GLpApp::GLpYUEpetraDataPool>
264Teuchos::RCP<const Epetra_MultiVector>
272Teuchos::RCP<const Epetra_Map>
278Teuchos::RCP<const Epetra_Map>
284Teuchos::RCP<const Epetra_Map>
287 TEUCHOS_TEST_FOR_EXCEPT(!(0<=l<=
Np_));
291Teuchos::RCP<const Epetra_Map>
294 TEUCHOS_TEST_FOR_EXCEPT(j!=0);
298Teuchos::RCP<const Epetra_Vector>
304Teuchos::RCP<const Epetra_Vector>
307 TEUCHOS_TEST_FOR_EXCEPT(!(0<=l<=
Np_));
311Teuchos::RCP<const Epetra_Vector>
317Teuchos::RCP<const Epetra_Vector>
323Teuchos::RCP<const Epetra_Vector>
326 TEUCHOS_TEST_FOR_EXCEPT(!(0<=l<=
Np_));
330Teuchos::RCP<const Epetra_Vector>
333 TEUCHOS_TEST_FOR_EXCEPT(!(0<=l<=
Np_));
337Teuchos::RCP<Epetra_Operator>
343Teuchos::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); } }
579 DfDp_op = &dyn_cast<Epetra_CrsMatrix>(*DfDp_out.
getLinearOp());
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_MultiVector > getMultiVector() const
Simple aggregate class that stores a derivative object as a general linear operator or as a multi-vec...
DerivativeMultiVector getDerivativeMultiVector() const
Teuchos::RCP< Epetra_Operator > getLinearOp() const
void setSupports(EInArgsMembers arg, bool supports=true)
void setModelEvalDescription(const std::string &modelEvalDescription)
Teuchos::RCP< const Epetra_Vector > get_p(int l) const
Teuchos::RCP< const Epetra_Vector > get_x() const
Set solution vector Taylor polynomial.
void setModelEvalDescription(const std::string &modelEvalDescription)
void set_W_properties(const DerivativeProperties &properties)
void set_DfDp_properties(int l, const DerivativeProperties &properties)
void set_DgDp_properties(int j, int l, const DerivativeProperties &properties)
void set_DgDx_properties(int j, const DerivativeProperties &properties)
void set_Np_Ng(int Np, int Ng)
void setSupports(EOutArgsMembers arg, bool supports=true)
Teuchos::RCP< Epetra_Operator > get_W() const
Derivative get_DfDp(int l) const
Evaluation< Epetra_Vector > get_g(int j) const
Get g(j) where 0 <= j && j < this->Ng().
Evaluation< Epetra_Vector > get_f() const
@ DERIV_LINEARITY_NONCONST
int ExtractMyRowView(int MyRow, int &NumEntries, double *&Values, int *&Indices) const
virtual void Print(std::ostream &os) const
int Dot(const Epetra_MultiVector &A, double *Result) const
virtual void Print(std::ostream &os) const
int Multiply(char TransA, char TransB, double ScalarAB, const Epetra_MultiVector &A, const Epetra_MultiVector &B, double ScalarThis)
int Update(double ScalarA, const Epetra_MultiVector &A, double ScalarThis)
int PutScalar(double ScalarConstant)
Teuchos::RCP< Epetra_Operator > create_W() const
Teuchos::RCP< Epetra_Vector > xU_
Teuchos::RCP< const Epetra_Map > get_p_map(int l) const
\breif .
InArgs createInArgs() const
Teuchos::RCP< Epetra_MultiVector > B_bar_
Teuchos::RCP< const Epetra_Map > map_x_
static const int p_rx_idx
RCP_Eptra_Map_Array_t map_p_
Teuchos::RCP< const Epetra_Map > map_p_bar_
Teuchos::RCP< const Epetra_Vector > get_x_init() const
Teuchos::RCP< const Epetra_Vector > get_p_upper_bounds(int l) const
Teuchos::RCP< Epetra_Vector > xL_
RCP_Eptra_Vector_Array_t p0_
Teuchos::RCP< GLpApp::GLpYUEpetraDataPool > dat_
Teuchos::RCP< const Epetra_Map > get_g_map(int j) const
\breif .
RCP_Eptra_Vector_Array_t pU_
Teuchos::RCP< const Epetra_Vector > q_
Teuchos::RCP< const Epetra_Map > map_f_
Teuchos::RCP< Epetra_CrsGraph > W_graph_
Teuchos::RCP< GLpApp::GLpYUEpetraDataPool > getDataPool()
Teuchos::RCP< const Epetra_Vector > get_x_upper_bounds() const
Teuchos::RCP< const Epetra_Vector > get_p_lower_bounds(int l) const
Teuchos::RCP< Epetra_Operator > create_DfDp_op(int l) const
Teuchos::RCP< const Epetra_MultiVector > get_B_bar() const
OutArgs createOutArgs() const
Teuchos::RCP< const Epetra_Vector > get_p_init(int l) const
Teuchos::RCP< const Epetra_Map > get_f_map() const
void set_q(Teuchos::RCP< const Epetra_Vector > const &q)
Teuchos::RCP< const Epetra_Map > map_g_
void evalModel(const InArgs &inArgs, const OutArgs &outArgs) 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< const Epetra_Vector > get_x_lower_bounds() const
static const int p_bndy_idx
Teuchos::RCP< const Epetra_Map > get_x_map() const
RCP_Eptra_Vector_Array_t pL_
Teuchos::RCP< Epetra_Vector > x0_