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 82 inline double sqr(
const double& s ) {
return s*s; }
102 ,
const char meshFile[]
106 ,
const double reactionRate
107 ,
const bool normalizeBasis
108 ,
const bool supportDerivatives
110 : supportDerivatives_(supportDerivatives), np_(np)
113 #ifdef GLPAPP_ADVDIFFREACT_OPTMODEL_DUMP_STUFF 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);
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 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 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";
241 #ifdef GLPAPP_ADVDIFFREACT_OPTMODEL_DUMP_STUFF 242 *out <<
"\nLeaving AdvDiffReactOptModel::AdvDiffReactOptModel(...) ...\n";
249 #ifdef GLPAPP_ADVDIFFREACT_OPTMODEL_DUMP_STUFF 253 *dout <<
"\nIn AdvDiffReactOptModel::set_q(q): q =\n\n";
414 #ifdef GLPAPP_ADVDIFFREACT_OPTMODEL_DUMP_STUFF 418 *dout <<
"\nEntering AdvDiffReactOptModel::evalModel(...) ...\n";
421 using Teuchos::rcp_dynamic_cast;
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";
440 *dout <<
"\np =\n\n";
474 if( g_out || DgDp_trans_out ) {
492 f.Update(1.0,Ax,0.0);
493 if(reactionRate!=0.0) {
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 517 *dout <<
"x-q =\n\n";
519 *dout <<
"H*(x-q) =\n\n";
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";
531 if(reactionRate!=0.0)
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;
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);
551 for(
int k = 0; k < DfDx_num_row_entries; ++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) {
563 DfDx_row_vals[k] = dat_A_row_vals[k];
569 if(out.
get() && trace) *out <<
"\nComputing DfDp ...\n";
576 if(out.
get() && dumpAll)
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);
615 for(
int k = 0; k < DfDp_num_row_entries; ++k) {
619 DfDp_row_vals[k] = dat_B_row_vals[k];
635 thyra_etaVec = Thyra::create_Vector(etaVec,space_p_bar);
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";
649 *dout <<
"\nDfDp_mv =\n\n";
661 xq.Update(-1.0,*
q_,1.0);
663 #ifdef GLPAPP_ADVDIFFREACT_OPTMODEL_DUMP_STUFF 668 *dout <<
"x-q =\n\n";
670 *dout <<
"DgDx_trans = H*(x-q) =\n\n";
686 #ifdef GLPAPP_ADVDIFFREACT_OPTMODEL_DUMP_STUFF 687 *dout <<
"\nR_p_bar =\n\n";
690 *dout <<
"\nB_bar =\n\n";
693 *dout <<
"\nDgDp_trans =\n\n";
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_FECrsMatrix > getR()
Teuchos::RCP< Epetra_FEVector > getq()
int MyGlobalElements(int *MyGlobalElementList) const
Teuchos::RCP< Epetra_Vector > x0_
Teuchos::RCP< const Epetra_IntSerialDenseVector > getpindx()
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)
const Epetra_CrsGraph & Graph() const
Teuchos::RCP< Epetra_Operator > create_DfDp_op(int l) const
Teuchos::RCP< const Epetra_Map > get_x_map() const
int Multiply(bool TransA, const Epetra_Vector &x, Epetra_Vector &y) const
void computeNpy(const Teuchos::RCP< const Epetra_MultiVector > &y)
Calls the function that computes the Jacobian of the nonlinear term.
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)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
virtual EVerbosityLevel getVerbLevel() const
void computeNy(const Teuchos::RCP< const Epetra_MultiVector > &y)
Calls the function that computes the nonlinear term.
Teuchos::RCP< const Epetra_Vector > get_p_lower_bounds(int l) const
void setSupports(EInArgsMembers arg, bool supports=true)
Teuchos::RCP< GLpApp::GLpYUEpetraDataPool > dat_
virtual RCP< FancyOStream > getOStream() const
InArgs createInArgs() const
Simple aggregate class that stores a derivative object as a general linear operator or as a multi-vec...
virtual std::string description() const
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_FECrsMatrix > getNpy()
Teuchos::RCP< Epetra_Operator > create_W() const
const Epetra_Map & OperatorDomainMap() 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
.
static RCP< Time > getNewTimer(const std::string &name)
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_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Teuchos::RCP< Epetra_Operator > get_W() const
Teuchos::RCP< const Epetra_Map > map_x_
T_To & dyn_cast(T_From &from)
int NumMyElements() const
static RCP< FancyOStream > getDefaultOStream()
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< Epetra_FECrsMatrix > getH()
Teuchos::RCP< const Epetra_Vector > get_p_init(int l) const
void setModelEvalDescription(const std::string &modelEvalDescription)
Teuchos::RCP< const Epetra_SerialDenseMatrix > getipcoords()
const Epetra_Map & OperatorRangeMap() const
Teuchos::RCP< Epetra_FECrsMatrix > getB()
void resize(size_type new_size, const value_type &x=value_type())
DerivativeMultiVector getDerivativeMultiVector() const
Teuchos::RCP< const Epetra_Vector > get_p(int l) const
virtual void Print(std::ostream &os) const
Teuchos::RCP< Epetra_FECrsMatrix > getA()
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
int NumGlobalElements() 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< Epetra_FEVector > getNy()
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)
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
Teuchos::RCP< Epetra_FEVector > getb()
Teuchos::RCP< const Epetra_Map > map_p_bar_