54 #include "ROL_TeuchosBatchManager.hpp"
56 #include "Teuchos_LAPACK.hpp"
84 void update(std::vector<Real> &u,
const std::vector<Real> &s,
const Real alpha=1.0)
const {
85 for (
unsigned i=0; i<u.size(); i++) {
90 void axpy(std::vector<Real> &out,
const Real a,
const std::vector<Real> &x,
const std::vector<Real> &y)
const {
91 for (
unsigned i=0; i < x.size(); i++) {
92 out[i] = a*x[i] + y[i];
96 void scale(std::vector<Real> &u,
const Real alpha=0.0)
const {
97 for (
unsigned i=0; i<u.size(); i++) {
102 void linear_solve(std::vector<Real> &u, std::vector<Real> &dl, std::vector<Real> &d, std::vector<Real> &du,
103 const std::vector<Real> &r,
const bool transpose =
false)
const {
104 if ( r.size() == 1 ) {
105 u.resize(1,r[0]/d[0]);
107 else if ( r.size() == 2 ) {
109 Real det = d[0]*d[1] - dl[0]*du[0];
110 u[0] = (d[1]*r[0] - du[0]*r[1])/det;
111 u[1] = (d[0]*r[1] - dl[0]*r[0])/det;
114 u.assign(r.begin(),r.end());
116 Teuchos::LAPACK<int,Real> lp;
117 std::vector<Real> du2(r.size()-2,0.0);
118 std::vector<int> ipiv(r.size(),0);
123 lp.GTTRF(dim,&dl[0],&d[0],&du[0],&du2[0],&ipiv[0],&info);
128 lp.GTTRS(trans,dim,nhrs,&dl[0],&d[0],&du[0],&du2[0],&ipiv[0],&u[0],ldb,&info);
133 BurgersFEM(
int nx = 128, Real nl = 1.0, Real cH1 = 1.0, Real cL2 = 1.0)
137 nu_ = std::pow(10.0,nu-2.0);
159 Real
compute_L2_dot(
const std::vector<Real> &x,
const std::vector<Real> &y)
const {
161 Real c = (((int)x.size()==
nx_) ? 4.0 : 2.0);
162 for (
unsigned i=0; i<x.size(); i++) {
164 ip +=
dx_/6.0*(c*x[i] + x[i+1])*y[i];
166 else if ( i == x.size()-1 ) {
167 ip +=
dx_/6.0*(x[i-1] + c*x[i])*y[i];
170 ip +=
dx_/6.0*(x[i-1] + 4.0*x[i] + x[i+1])*y[i];
182 void apply_mass(std::vector<Real> &Mu,
const std::vector<Real> &u )
const {
183 Mu.resize(u.size(),0.0);
184 Real c = (((int)u.size()==
nx_) ? 4.0 : 2.0);
185 for (
unsigned i=0; i<u.size(); i++) {
187 Mu[i] =
dx_/6.0*(c*u[i] + u[i+1]);
189 else if ( i == u.size()-1 ) {
190 Mu[i] =
dx_/6.0*(u[i-1] + c*u[i]);
193 Mu[i] =
dx_/6.0*(u[i-1] + 4.0*u[i] + u[i+1]);
200 unsigned nx = u.size();
202 std::vector<Real> dl(nx-1,
dx_/6.0);
203 std::vector<Real> d(nx,2.0*
dx_/3.0);
204 std::vector<Real> du(nx-1,
dx_/6.0);
205 if ( (
int)nx !=
nx_ ) {
214 std::vector<Real> u(
nx_,0.0), Mu(
nx_,0.0), iMMu(
nx_,0.0), diff(
nx_,0.0);
215 for (
int i = 0; i <
nx_; i++) {
216 u[i] = 2.0*(Real)rand()/(Real)RAND_MAX - 1.0;
220 axpy(diff,-1.0,iMMu,u);
223 outStream <<
"Test Inverse State Mass Matrix\n";
224 outStream <<
" ||u - inv(M)Mu|| = " << error <<
"\n";
225 outStream <<
" ||u|| = " << normu <<
"\n";
226 outStream <<
" Relative Error = " << error/normu <<
"\n";
229 u.resize(nx_+2,0.0); Mu.resize(nx_+2,0.0); iMMu.resize(nx_+2,0.0); diff.resize(nx_+2,0.0);
230 for (
int i = 0; i < nx_+2; i++) {
231 u[i] = 2.0*(Real)rand()/(Real)RAND_MAX - 1.0;
235 axpy(diff,-1.0,iMMu,u);
238 outStream <<
"Test Inverse Control Mass Matrix\n";
239 outStream <<
" ||z - inv(M)Mz|| = " << error <<
"\n";
240 outStream <<
" ||z|| = " << normu <<
"\n";
241 outStream <<
" Relative Error = " << error/normu <<
"\n";
249 Real
compute_H1_dot(
const std::vector<Real> &x,
const std::vector<Real> &y)
const {
251 for (
int i=0; i<
nx_; i++) {
253 ip +=
cL2_*
dx_/6.0*(4.0*x[i] + x[i+1])*y[i];
254 ip +=
cH1_*(2.0*x[i] - x[i+1])/
dx_*y[i];
256 else if ( i == nx_-1 ) {
257 ip +=
cL2_*
dx_/6.0*(x[i-1] + 4.0*x[i])*y[i];
258 ip +=
cH1_*(2.0*x[i] - x[i-1])/
dx_*y[i];
261 ip +=
cL2_*
dx_/6.0*(x[i-1] + 4.0*x[i] + x[i+1])*y[i];
262 ip +=
cH1_*(2.0*x[i] - x[i-1] - x[i+1])/
dx_*y[i];
274 void apply_H1(std::vector<Real> &Mu,
const std::vector<Real> &u )
const {
276 for (
int i=0; i<
nx_; i++) {
278 Mu[i] =
cL2_*
dx_/6.0*(4.0*u[i] + u[i+1])
279 +
cH1_*(2.0*u[i] - u[i+1])/
dx_;
281 else if ( i == nx_-1 ) {
282 Mu[i] =
cL2_*
dx_/6.0*(u[i-1] + 4.0*u[i])
283 +
cH1_*(2.0*u[i] - u[i-1])/
dx_;
286 Mu[i] =
cL2_*
dx_/6.0*(u[i-1] + 4.0*u[i] + u[i+1])
287 +
cH1_*(2.0*u[i] - u[i-1] - u[i+1])/
dx_;
302 std::vector<Real> u(
nx_,0.0), Mu(
nx_,0.0), iMMu(
nx_,0.0), diff(
nx_,0.0);
303 for (
int i = 0; i <
nx_; i++) {
304 u[i] = 2.0*(Real)rand()/(Real)RAND_MAX - 1.0;
308 axpy(diff,-1.0,iMMu,u);
311 outStream <<
"Test Inverse State H1 Matrix\n";
312 outStream <<
" ||u - inv(M)Mu|| = " << error <<
"\n";
313 outStream <<
" ||u|| = " << normu <<
"\n";
314 outStream <<
" Relative Error = " << error/normu <<
"\n";
323 const std::vector<Real> &z)
const {
326 for (
int i=0; i<
nx_; i++) {
329 r[i] =
nu_/
dx_*(2.0*u[i]-u[i+1]);
332 r[i] =
nu_/
dx_*(2.0*u[i]-u[i-1]);
335 r[i] =
nu_/
dx_*(2.0*u[i]-u[i-1]-u[i+1]);
339 r[i] +=
nl_*u[i+1]*(u[i]+u[i+1])/6.0;
342 r[i] -=
nl_*u[i-1]*(u[i-1]+u[i])/6.0;
345 r[i] -=
dx_/6.0*(z[i]+4.0*z[i+1]+z[i+2]);
359 std::vector<Real> &d,
360 std::vector<Real> &du,
361 const std::vector<Real> &u)
const {
370 for (
int i=0; i<
nx_; i++) {
372 dl[i] +=
nl_*(-2.0*u[i]-u[i+1])/6.0;
373 d[i] +=
nl_*u[i+1]/6.0;
376 d[i] -=
nl_*u[i-1]/6.0;
377 du[i-1] +=
nl_*(u[i-1]+2.0*u[i])/6.0;
387 const std::vector<Real> &v,
388 const std::vector<Real> &u,
389 const std::vector<Real> &z)
const {
391 for (
int i = 0; i <
nx_; i++) {
394 jv[i] += -
nu_/
dx_*v[i-1]-
nl_*(u[i-1]/6.0*v[i]+(u[i]+2.0*u[i-1])/6.0*v[i-1]);
397 jv[i] += -
nu_/
dx_*v[i+1]+
nl_*(u[i+1]/6.0*v[i]+(u[i]+2.0*u[i+1])/6.0*v[i+1]);
400 jv[ 0] -=
nl_*
u0_/6.0*v[0];
401 jv[nx_-1] +=
nl_*
u1_/6.0*v[nx_-1];
406 const std::vector<Real> &v,
407 const std::vector<Real> &u,
408 const std::vector<Real> &z)
const {
410 std::vector<Real> d(
nx_,0.0);
411 std::vector<Real> dl(
nx_-1,0.0);
412 std::vector<Real> du(
nx_-1,0.0);
420 const std::vector<Real> &v,
421 const std::vector<Real> &u,
422 const std::vector<Real> &z)
const {
424 for (
int i = 0; i <
nx_; i++) {
425 ajv[i] =
nu_/
dx_*2.0*v[i];
427 ajv[i] += -
nu_/
dx_*v[i-1]-
nl_*(u[i-1]/6.0*v[i]
428 -(u[i-1]+2.0*u[i])/6.0*v[i-1]);
431 ajv[i] += -
nu_/
dx_*v[i+1]+
nl_*(u[i+1]/6.0*v[i]
432 -(u[i+1]+2.0*u[i])/6.0*v[i+1]);
435 ajv[ 0] -=
nl_*
u0_/6.0*v[0];
436 ajv[nx_-1] +=
nl_*
u1_/6.0*v[nx_-1];
441 const std::vector<Real> &v,
442 const std::vector<Real> &u,
443 const std::vector<Real> &z)
const {
445 std::vector<Real> d(
nx_,0.0);
446 std::vector<Real> du(
nx_-1,0.0);
447 std::vector<Real> dl(
nx_-1,0.0);
458 const std::vector<Real> &v,
459 const std::vector<Real> &u,
460 const std::vector<Real> &z)
const {
461 for (
int i=0; i<
nx_; i++) {
463 jv[i] = -
dx_/6.0*(v[i]+4.0*v[i+1]+v[i+2]);
469 const std::vector<Real> &v,
470 const std::vector<Real> &u,
471 const std::vector<Real> &z)
const {
472 for (
int i=0; i<
nx_+2; i++) {
474 jv[i] = -
dx_/6.0*v[i];
477 jv[i] = -
dx_/6.0*(4.0*v[i-1]+v[i]);
479 else if ( i == nx_ ) {
480 jv[i] = -
dx_/6.0*(4.0*v[i-1]+v[i-2]);
482 else if ( i == nx_+1 ) {
483 jv[i] = -
dx_/6.0*v[i-2];
486 jv[i] = -
dx_/6.0*(v[i-2]+4.0*v[i-1]+v[i]);
495 const std::vector<Real> &w,
496 const std::vector<Real> &v,
497 const std::vector<Real> &u,
498 const std::vector<Real> &z)
const {
499 for (
int i=0; i<
nx_; i++) {
503 ahwv[i] += (w[i]*v[i+1] - w[i+1]*(2.0*v[i]+v[i+1]))/6.0;
506 ahwv[i] += (w[i-1]*(v[i-1]+2.0*v[i]) - w[i]*v[i-1])/6.0;
512 const std::vector<Real> &w,
513 const std::vector<Real> &v,
514 const std::vector<Real> &u,
515 const std::vector<Real> &z) {
516 ahwv.assign(u.size(),0.0);
519 const std::vector<Real> &w,
520 const std::vector<Real> &v,
521 const std::vector<Real> &u,
522 const std::vector<Real> &z) {
523 ahwv.assign(z.size(),0.0);
526 const std::vector<Real> &w,
527 const std::vector<Real> &v,
528 const std::vector<Real> &u,
529 const std::vector<Real> &z) {
530 ahwv.assign(z.size(),0.0);
537 ROL::Ptr<std::vector<Real> >
vec_;
538 ROL::Ptr<BurgersFEM<Real> >
fem_;
540 mutable ROL::Ptr<L2VectorDual<Real> >
dual_vec_;
549 const std::vector<Real>& xval = *ex.
getVector();
550 std::copy(xval.begin(),xval.end(),
vec_->begin());
555 const std::vector<Real>& xval = *ex.
getVector();
558 (*vec_)[i] += xval[i];
571 const std::vector<Real>& xval = *ex.
getVector();
572 return fem_->compute_L2_dot(xval,*
vec_);
577 val = std::sqrt(
dot(*
this) );
581 ROL::Ptr<ROL::Vector<Real> >
clone()
const {
582 return ROL::makePtr<L2VectorPrimal>( ROL::makePtr<std::vector<Real>>(
vec_->size(),0.0),
fem_);
593 ROL::Ptr<ROL::Vector<Real> >
basis(
const int i )
const {
594 ROL::Ptr<L2VectorPrimal> e
595 = ROL::makePtr<L2VectorPrimal>( ROL::makePtr<std::vector<Real>>(
vec_->size(),0.0),
fem_);
596 (*e->getVector())[i] = 1.0;
605 dual_vec_ = ROL::makePtr<L2VectorDual<Real>>(
606 ROL::makePtr<std::vector<Real>>(*vec_),
fem_);
608 fem_->apply_mass(*(ROL::constPtrCast<std::vector<Real> >(
dual_vec_->getVector())),*
vec_);
614 const std::vector<Real>& xval = *ex.
getVector();
615 return std::inner_product(
vec_->begin(),
vec_->end(), xval.begin(), Real(0));
623 ROL::Ptr<std::vector<Real> >
vec_;
624 ROL::Ptr<BurgersFEM<Real> >
fem_;
626 mutable ROL::Ptr<L2VectorPrimal<Real> >
dual_vec_;
635 const std::vector<Real>& xval = *ex.
getVector();
636 std::copy(xval.begin(),xval.end(),
vec_->begin());
641 const std::vector<Real>& xval = *ex.
getVector();
644 (*vec_)[i] += xval[i];
657 const std::vector<Real>& xval = *ex.
getVector();
659 std::vector<Real> Mx(dimension,0.0);
660 fem_->apply_inverse_mass(Mx,xval);
662 for (
unsigned i = 0; i <
dimension; i++) {
663 val += Mx[i]*(*vec_)[i];
670 val = std::sqrt(
dot(*
this) );
674 ROL::Ptr<ROL::Vector<Real> >
clone()
const {
675 return ROL::makePtr<L2VectorDual>( ROL::makePtr<std::vector<Real>>(
vec_->size(),0.0),
fem_);
686 ROL::Ptr<ROL::Vector<Real> >
basis(
const int i )
const {
687 ROL::Ptr<L2VectorDual> e
688 = ROL::makePtr<L2VectorDual>( ROL::makePtr<std::vector<Real>>(
vec_->size(),0.0),
fem_);
689 (*e->getVector())[i] = 1.0;
698 dual_vec_ = ROL::makePtr<L2VectorPrimal<Real>>(
699 ROL::makePtr<std::vector<Real>>(*vec_),
fem_);
701 fem_->apply_inverse_mass(*(ROL::constPtrCast<std::vector<Real> >(
dual_vec_->getVector())),*
vec_);
707 const std::vector<Real>& xval = *ex.
getVector();
708 return std::inner_product(
vec_->begin(),
vec_->end(), xval.begin(), Real(0));
716 ROL::Ptr<std::vector<Real> >
vec_;
717 ROL::Ptr<BurgersFEM<Real> >
fem_;
719 mutable ROL::Ptr<H1VectorDual<Real> >
dual_vec_;
728 const std::vector<Real>& xval = *ex.
getVector();
729 std::copy(xval.begin(),xval.end(),
vec_->begin());
734 const std::vector<Real>& xval = *ex.
getVector();
737 (*vec_)[i] += xval[i];
750 const std::vector<Real>& xval = *ex.
getVector();
751 return fem_->compute_H1_dot(xval,*
vec_);
756 val = std::sqrt(
dot(*
this) );
760 ROL::Ptr<ROL::Vector<Real> >
clone()
const {
761 return ROL::makePtr<H1VectorPrimal>( ROL::makePtr<std::vector<Real>>(
vec_->size(),0.0),
fem_);
772 ROL::Ptr<ROL::Vector<Real> >
basis(
const int i )
const {
773 ROL::Ptr<H1VectorPrimal> e
774 = ROL::makePtr<H1VectorPrimal>( ROL::makePtr<std::vector<Real>>(
vec_->size(),0.0),
fem_);
775 (*e->getVector())[i] = 1.0;
784 dual_vec_ = ROL::makePtr<H1VectorDual<Real>>(
785 ROL::makePtr<std::vector<Real>>(*vec_),
fem_);
787 fem_->apply_H1(*(ROL::constPtrCast<std::vector<Real> >(
dual_vec_->getVector())),*
vec_);
793 const std::vector<Real>& xval = *ex.
getVector();
794 return std::inner_product(
vec_->begin(),
vec_->end(), xval.begin(), Real(0));
802 ROL::Ptr<std::vector<Real> >
vec_;
803 ROL::Ptr<BurgersFEM<Real> >
fem_;
805 mutable ROL::Ptr<H1VectorPrimal<Real> >
dual_vec_;
814 const std::vector<Real>& xval = *ex.
getVector();
815 std::copy(xval.begin(),xval.end(),
vec_->begin());
820 const std::vector<Real>& xval = *ex.
getVector();
823 (*vec_)[i] += xval[i];
836 const std::vector<Real>& xval = *ex.
getVector();
838 std::vector<Real> Mx(dimension,0.0);
839 fem_->apply_inverse_H1(Mx,xval);
841 for (
unsigned i = 0; i <
dimension; i++) {
842 val += Mx[i]*(*vec_)[i];
849 val = std::sqrt(
dot(*
this) );
853 ROL::Ptr<ROL::Vector<Real> >
clone()
const {
854 return ROL::makePtr<H1VectorDual>( ROL::makePtr<std::vector<Real>>(
vec_->size(),0.0),
fem_);
865 ROL::Ptr<ROL::Vector<Real> >
basis(
const int i )
const {
866 ROL::Ptr<H1VectorDual> e
867 = ROL::makePtr<H1VectorDual>( ROL::makePtr<std::vector<Real>>(
vec_->size(),0.0),
fem_);
868 (*e->getVector())[i] = 1.0;
877 dual_vec_ = ROL::makePtr<H1VectorPrimal<Real>>(
878 ROL::makePtr<std::vector<Real>>(*vec_),
fem_);
880 fem_->apply_inverse_H1(*(ROL::constPtrCast<std::vector<Real> >(
dual_vec_->getVector())),*
vec_);
886 const std::vector<Real>& xval = *ex.
getVector();
887 return std::inner_product(
vec_->begin(),
vec_->end(), xval.begin(), Real(0));
907 ROL::Ptr<BurgersFEM<Real> >
fem_;
916 ROL::Ptr<std::vector<Real> > cp =
918 ROL::Ptr<const std::vector<Real> > up =
920 ROL::Ptr<const std::vector<Real> > zp =
923 const std::vector<Real> param
925 fem_->set_problem_data(param[0],param[1],param[2],param[3]);
927 fem_->compute_residual(*cp,*up,*zp);
932 ROL::Ptr<std::vector<Real> > jvp =
934 ROL::Ptr<const std::vector<Real> > vp =
936 ROL::Ptr<const std::vector<Real> > up =
938 ROL::Ptr<const std::vector<Real> > zp =
941 const std::vector<Real> param
943 fem_->set_problem_data(param[0],param[1],param[2],param[3]);
945 fem_->apply_pde_jacobian(*jvp,*vp,*up,*zp);
950 ROL::Ptr<std::vector<Real> > jvp =
952 ROL::Ptr<const std::vector<Real> > vp =
954 ROL::Ptr<const std::vector<Real> > up =
956 ROL::Ptr<const std::vector<Real> > zp =
959 const std::vector<Real> param
961 fem_->set_problem_data(param[0],param[1],param[2],param[3]);
963 fem_->apply_control_jacobian(*jvp,*vp,*up,*zp);
968 ROL::Ptr<std::vector<Real> > ijvp =
970 ROL::Ptr<const std::vector<Real> > vp =
972 ROL::Ptr<const std::vector<Real> > up =
974 ROL::Ptr<const std::vector<Real> > zp =
977 const std::vector<Real> param
979 fem_->set_problem_data(param[0],param[1],param[2],param[3]);
981 fem_->apply_inverse_pde_jacobian(*ijvp,*vp,*up,*zp);
986 ROL::Ptr<std::vector<Real> > jvp =
988 ROL::Ptr<const std::vector<Real> > vp =
990 ROL::Ptr<const std::vector<Real> > up =
992 ROL::Ptr<const std::vector<Real> > zp =
995 const std::vector<Real> param
997 fem_->set_problem_data(param[0],param[1],param[2],param[3]);
999 fem_->apply_adjoint_pde_jacobian(*jvp,*vp,*up,*zp);
1004 ROL::Ptr<std::vector<Real> > jvp =
1006 ROL::Ptr<const std::vector<Real> > vp =
1008 ROL::Ptr<const std::vector<Real> > up =
1010 ROL::Ptr<const std::vector<Real> > zp =
1013 const std::vector<Real> param
1015 fem_->set_problem_data(param[0],param[1],param[2],param[3]);
1017 fem_->apply_adjoint_control_jacobian(*jvp,*vp,*up,*zp);
1022 ROL::Ptr<std::vector<Real> > iajvp =
1024 ROL::Ptr<const std::vector<Real> > vp =
1026 ROL::Ptr<const std::vector<Real> > up =
1028 ROL::Ptr<const std::vector<Real> > zp =
1031 const std::vector<Real> param
1033 fem_->set_problem_data(param[0],param[1],param[2],param[3]);
1035 fem_->apply_inverse_adjoint_pde_jacobian(*iajvp,*vp,*up,*zp);
1041 ROL::Ptr<std::vector<Real> > ahwvp =
1043 ROL::Ptr<const std::vector<Real> > wp =
1045 ROL::Ptr<const std::vector<Real> > vp =
1047 ROL::Ptr<const std::vector<Real> > up =
1049 ROL::Ptr<const std::vector<Real> > zp =
1052 const std::vector<Real> param
1054 fem_->set_problem_data(param[0],param[1],param[2],param[3]);
1056 fem_->apply_adjoint_pde_hessian(*ahwvp,*wp,*vp,*up,*zp);
1066 ROL::Ptr<std::vector<Real> > ahwvp =
1068 ROL::Ptr<const std::vector<Real> > wp =
1070 ROL::Ptr<const std::vector<Real> > vp =
1072 ROL::Ptr<const std::vector<Real> > up =
1074 ROL::Ptr<const std::vector<Real> > zp =
1077 const std::vector<Real> param
1079 fem_->set_problem_data(param[0],param[1],param[2],param[3]);
1081 fem_->apply_adjoint_control_pde_hessian(*ahwvp,*wp,*vp,*up,*zp);
1090 ROL::Ptr<std::vector<Real> > ahwvp =
1092 ROL::Ptr<const std::vector<Real> > wp =
1094 ROL::Ptr<const std::vector<Real> > vp =
1096 ROL::Ptr<const std::vector<Real> > up =
1098 ROL::Ptr<const std::vector<Real> > zp =
1101 const std::vector<Real> param
1103 fem_->set_problem_data(param[0],param[1],param[2],param[3]);
1105 fem_->apply_adjoint_pde_control_hessian(*ahwvp,*wp,*vp,*up,*zp);
1114 ROL::Ptr<std::vector<Real> > ahwvp =
1116 ROL::Ptr<const std::vector<Real> > wp =
1118 ROL::Ptr<const std::vector<Real> > vp =
1120 ROL::Ptr<const std::vector<Real> > up =
1122 ROL::Ptr<const std::vector<Real> > zp =
1125 const std::vector<Real> param
1127 fem_->set_problem_data(param[0],param[1],param[2],param[3]);
1129 fem_->apply_adjoint_control_hessian(*ahwvp,*wp,*vp,*up,*zp);
1137 template<
class Real>
1149 ROL::Ptr<BurgersFEM<Real> >
fem_;
1156 Real x = 0.0) :
fem_(fem),
x_(x) {
1157 for (
int i = 1; i <
fem_->num_dof()+1; i++) {
1158 if ( (Real)i*(
fem_->mesh_spacing()) >=
x_ ) {
1167 ROL::Ptr<const std::vector<Real> > up =
1177 Real val = 0.5*((((Real)
indices_[0]+1.)*(
fem_->mesh_spacing())-
x_)
1180 for (
uint i = 1; i < indices_.size(); i++) {
1181 val += (
fem_->mesh_spacing())*(*up)[indices_[i]];
1187 ROL::Ptr<std::vector<Real> > gp =
1189 ROL::Ptr<const std::vector<Real> > up =
1201 (*gp)[
indices_[0]] = -0.5*((((Real)indices_[0]+1.)*(
fem_->mesh_spacing())-
x_)
1202 *(
x_+(2.-((Real)indices_[0]+1.))*(
fem_->mesh_spacing()))/(
fem_->mesh_spacing())
1203 +(
fem_->mesh_spacing()));
1206 for (
uint i = 1; i < indices_.size(); i++) {
1207 (*gp)[indices_[i]] = -(
fem_->mesh_spacing());
1249 template<
class Real>
1253 std::vector<Real>
x_lo_;
1254 std::vector<Real>
x_up_;
1257 ROL::Ptr<BurgersFEM<Real> >
fem_;
1258 ROL::Ptr<ROL::Vector<Real> >
l_;
1259 ROL::Ptr<ROL::Vector<Real> >
u_;
1266 catch (std::exception &e) {
1276 catch (std::exception &e) {
1281 void axpy(std::vector<Real> &out,
const Real a,
1282 const std::vector<Real> &x,
const std::vector<Real> &y)
const{
1283 out.resize(
dim_,0.0);
1284 for (
unsigned i = 0; i <
dim_; i++) {
1285 out[i] = a*x[i] + y[i];
1290 for (
int i = 0; i <
dim_; i++ ) {
1291 x[i] = std::max(
x_lo_[i],std::min(
x_up_[i],x[i]));
1300 for (
int i = 0; i <
dim_; i++ ) {
1309 l_ = ROL::makePtr<L2VectorPrimal<Real>>(
1310 ROL::makePtr<std::vector<Real>>(l), fem);
1311 u_ = ROL::makePtr<L2VectorPrimal<Real>>(
1312 ROL::makePtr<std::vector<Real>>(u), fem);
1319 for (
int i = 0; i <
dim_; i++ ) {
1320 if ( (*ex)[i] >=
x_lo_[i] && (*ex)[i] <=
x_up_[i] ) { cnt *= 1; }
1323 if ( cnt == 0 ) { val =
false; }
1328 ROL::Ptr<std::vector<Real> > ex;
cast_vector(ex,x);
1334 ROL::Ptr<std::vector<Real> > ev;
cast_vector(ev,v);
1336 for (
int i = 0; i <
dim_; i++ ) {
1337 if ( ((*ex)[i] <=
x_lo_[i]+epsn) ) {
1345 ROL::Ptr<std::vector<Real> > ev;
cast_vector(ev,v);
1347 for (
int i = 0; i <
dim_; i++ ) {
1348 if ( ((*ex)[i] >=
x_up_[i]-epsn) ) {
1356 ROL::Ptr<std::vector<Real> > ev;
cast_vector(ev,v);
1358 for (
int i = 0; i <
dim_; i++ ) {
1359 if ( ((*ex)[i] <=
x_lo_[i]+epsn) ||
1360 ((*ex)[i] >=
x_up_[i]-epsn) ) {
1369 ROL::Ptr<std::vector<Real> > ev;
cast_vector(ev,v);
1371 for (
int i = 0; i <
dim_; i++ ) {
1372 if ( ((*ex)[i] <=
x_lo_[i]+epsn && (*eg)[i] > geps) ) {
1381 ROL::Ptr<std::vector<Real> > ev;
cast_vector(ev,v);
1383 for (
int i = 0; i <
dim_; i++ ) {
1384 if ( ((*ex)[i] >=
x_up_[i]-epsn && (*eg)[i] < -geps) ) {
1393 ROL::Ptr<std::vector<Real> > ev;
cast_vector(ev,v);
1395 for (
int i = 0; i <
dim_; i++ ) {
1396 if ( ((*ex)[i] <=
x_lo_[i]+epsn && (*eg)[i] > geps) ||
1397 ((*ex)[i] >=
x_up_[i]-epsn && (*eg)[i] < -geps) ) {
1412 template<
class Real>
1416 std::vector<Real>
x_lo_;
1417 std::vector<Real>
x_up_;
1420 ROL::Ptr<BurgersFEM<Real> >
fem_;
1421 ROL::Ptr<ROL::Vector<Real> >
l_;
1422 ROL::Ptr<ROL::Vector<Real> >
u_;
1427 xvec = ROL::constPtrCast<std::vector<Real> >(
1430 catch (std::exception &e) {
1431 xvec = ROL::constPtrCast<std::vector<Real> >(
1442 catch (std::exception &e) {
1448 void axpy(std::vector<Real> &out,
const Real a,
1449 const std::vector<Real> &x,
const std::vector<Real> &y)
const{
1450 out.resize(
dim_,0.0);
1451 for (
unsigned i = 0; i <
dim_; i++) {
1452 out[i] = a*x[i] + y[i];
1457 for (
int i = 0; i <
dim_; i++ ) {
1458 x[i] = std::max(
x_lo_[i],std::min(
x_up_[i],x[i]));
1467 for (
int i = 0; i <
dim_; i++ ) {
1476 l_ = ROL::makePtr<H1VectorPrimal<Real>>(
1477 ROL::makePtr<std::vector<Real>>(l), fem);
1478 u_ = ROL::makePtr<H1VectorPrimal<Real>>(
1479 ROL::makePtr<std::vector<Real>>(u), fem);
1486 for (
int i = 0; i <
dim_; i++ ) {
1487 if ( (*ex)[i] >=
x_lo_[i] && (*ex)[i] <=
x_up_[i] ) { cnt *= 1; }
1490 if ( cnt == 0 ) { val =
false; }
1495 ROL::Ptr<std::vector<Real> > ex;
cast_vector(ex,x);
1501 ROL::Ptr<std::vector<Real> > ev;
cast_vector(ev,v);
1503 for (
int i = 0; i <
dim_; i++ ) {
1504 if ( ((*ex)[i] <=
x_lo_[i]+epsn) ) {
1512 ROL::Ptr<std::vector<Real> > ev;
cast_vector(ev,v);
1514 for (
int i = 0; i <
dim_; i++ ) {
1515 if ( ((*ex)[i] >=
x_up_[i]-epsn) ) {
1523 ROL::Ptr<std::vector<Real> > ev;
cast_vector(ev,v);
1525 for (
int i = 0; i <
dim_; i++ ) {
1526 if ( ((*ex)[i] <=
x_lo_[i]+epsn) ||
1527 ((*ex)[i] >=
x_up_[i]-epsn) ) {
1536 ROL::Ptr<std::vector<Real> > ev;
cast_vector(ev,v);
1538 for (
int i = 0; i <
dim_; i++ ) {
1539 if ( ((*ex)[i] <=
x_lo_[i]+epsn && (*eg)[i] > geps) ) {
1548 ROL::Ptr<std::vector<Real> > ev;
cast_vector(ev,v);
1550 for (
int i = 0; i <
dim_; i++ ) {
1551 if ( ((*ex)[i] >=
x_up_[i]-epsn && (*eg)[i] < -geps) ) {
1560 ROL::Ptr<std::vector<Real> > ev;
cast_vector(ev,v);
1562 for (
int i = 0; i <
dim_; i++ ) {
1563 if ( ((*ex)[i] <=
x_lo_[i]+epsn && (*eg)[i] > geps) ||
1564 ((*ex)[i] >=
x_up_[i]-epsn && (*eg)[i] < -geps) ) {
1579 template<
class Real,
class Ordinal>
1585 xvec = ROL::constPtrCast<std::vector<Real> >(
1588 catch (std::exception &e) {
1589 xvec = ROL::constPtrCast<std::vector<Real> >(
1596 : ROL::TeuchosBatchManager<Real,Ordinal>(comm) {}
1598 ROL::Ptr<std::vector<Real> > input_ptr;
1600 int dim_i = input_ptr->size();
1601 ROL::Ptr<std::vector<Real> > output_ptr;
1603 int dim_o = output_ptr->size();
1604 if ( dim_i != dim_o ) {
1605 std::cout <<
"L2VectorBatchManager: DIMENSION MISMATCH ON RANK "
1606 << ROL::TeuchosBatchManager<Real,Ordinal>::batchID() <<
"\n";
1609 ROL::TeuchosBatchManager<Real,Ordinal>::sumAll(&(*input_ptr)[0],&(*output_ptr)[0],dim_i);
1614 template<
class Real,
class Ordinal>
1620 xvec = ROL::constPtrCast<std::vector<Real> >(
1623 catch (std::exception &e) {
1624 xvec = ROL::constPtrCast<std::vector<Real> >(
1631 : ROL::TeuchosBatchManager<Real,Ordinal>(comm) {}
1633 ROL::Ptr<std::vector<Real> > input_ptr;
1635 int dim_i = input_ptr->size();
1636 ROL::Ptr<std::vector<Real> > output_ptr;
1638 int dim_o = output_ptr->size();
1639 if ( dim_i != dim_o ) {
1640 std::cout <<
"H1VectorBatchManager: DIMENSION MISMATCH ON RANK "
1641 << ROL::TeuchosBatchManager<Real,Ordinal>::batchID() <<
"\n";
1644 ROL::TeuchosBatchManager<Real,Ordinal>::sumAll(&(*input_ptr)[0],&(*output_ptr)[0],dim_i);
1649 template<
class Real>
1650 Real
random(
const ROL::Ptr<
const Teuchos::Comm<int> > &comm) {
1652 if ( Teuchos::rank<int>(*comm)==0 ) {
1653 val = (Real)rand()/(Real)RAND_MAX;
1655 Teuchos::broadcast<int,Real>(*comm,0,1,&val);
H1VectorPrimal< Real > DualConstraintVector
BurgersFEM(int nx=128, Real nl=1.0, Real cH1=1.0, Real cL2=1.0)
L2VectorDual(const ROL::Ptr< std::vector< Real > > &vec, const ROL::Ptr< BurgersFEM< Real > > &fem)
ROL::Ptr< std::vector< Real > > vec_
Provides the interface to evaluate simulation-based objective functions.
std::vector< Real > x_up_
typename PV< Real >::size_type size_type
Real compute_H1_dot(const std::vector< Real > &x, const std::vector< Real > &y) const
void applyAdjointJacobian_1(ROL::Vector< Real > &ajv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Apply the adjoint of the partial constraint Jacobian at , , to the vector . This is the primary inter...
Real norm() const
Returns where .
void apply_mass(std::vector< Real > &Mu, const std::vector< Real > &u) const
bool isFeasible(const ROL::Vector< Real > &x)
Check if the vector, v, is feasible.
void pruneActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &g, const ROL::Vector< Real > &x, Real xeps=Real(0), Real geps=Real(0))
void axpy(std::vector< Real > &out, const Real a, const std::vector< Real > &x, const std::vector< Real > &y) const
ROL::Ptr< ROL::Vector< Real > > basis(const int i) const
Return i-th basis vector.
void cast_vector(ROL::Ptr< std::vector< Real > > &xvec, ROL::Vector< Real > &x) const
Real dot(const ROL::Vector< Real > &x) const override
Compute where .
Real dot(const ROL::Vector< Real > &x) const
Compute where .
void apply_adjoint_pde_control_hessian(std::vector< Real > &ahwv, const std::vector< Real > &w, const std::vector< Real > &v, const std::vector< Real > &u, const std::vector< Real > &z)
ROL::Ptr< std::vector< Real > > vec_
int dimension() const
Return dimension of the vector space.
ROL::Ptr< L2VectorDual< Real > > dual_vec_
Real get_viscosity(void) const
void applyAdjointHessian_12(ROL::Vector< Real > &ahwv, const ROL::Vector< Real > &w, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Apply the optimization-space derivative of the adjoint of the constraint simulation-space Jacobian at...
void plus(const ROL::Vector< Real > &x)
Compute , where .
std::vector< Real > x_up_
void pruneUpperActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &g, const ROL::Vector< Real > &x, Real xeps=Real(0), Real geps=Real(0))
Set variables to zero if they correspond to the upper -binding set.
void apply_adjoint_pde_hessian(std::vector< Real > &ahwv, const std::vector< Real > &w, const std::vector< Real > &v, const std::vector< Real > &u, const std::vector< Real > &z) const
ROL::Ptr< BurgersFEM< Real > > fem_
H1VectorPrimal(const ROL::Ptr< std::vector< Real > > &vec, const ROL::Ptr< BurgersFEM< Real > > &fem)
Real norm() const
Returns where .
Contains definitions of custom data types in ROL.
Real dot(const ROL::Vector< Real > &x) const
Compute where .
void cast_vector(ROL::Ptr< std::vector< Real > > &xvec, ROL::Vector< Real > &x) const
void applyAdjointHessian_11(ROL::Vector< Real > &ahwv, const ROL::Vector< Real > &w, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Apply the simulation-space derivative of the adjoint of the constraint simulation-space Jacobian at ...
ROL::Ptr< const std::vector< Real > > getVector() const
const ROL::Ptr< const ROL::Vector< Real > > getUpperBound(void) const
Return the ref count pointer to the upper bound vector.
void plus(const ROL::Vector< Real > &x)
Compute , where .
ROL::Ptr< BurgersFEM< Real > > fem_
const std::vector< Real > getParameter(void) const
Real compute_H1_norm(const std::vector< Real > &r) const
L2BoundConstraint(std::vector< Real > &l, std::vector< Real > &u, const ROL::Ptr< BurgersFEM< Real > > &fem, Real scale=1.0)
Real norm() const
Returns where .
Real apply(const ROL::Vector< Real > &x) const
Apply to a dual vector. This is equivalent to the call .
void applyInverseAdjointJacobian_1(ROL::Vector< Real > &iajv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Apply the inverse of the adjoint of the partial constraint Jacobian at , , to the vector ...
ROL::Ptr< std::vector< Real > > vec_
void pruneLowerActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &x, Real eps=Real(0))
Set variables to zero if they correspond to the lower -active set.
H1VectorDual(const ROL::Ptr< std::vector< Real > > &vec, const ROL::Ptr< BurgersFEM< Real > > &fem)
void project(ROL::Vector< Real > &x)
Project optimization variables onto the bounds.
virtual void zero()
Set to zero vector.
ROL::Ptr< std::vector< Real > > getVector()
Defines the linear algebra or vector space interface.
void cast_vector(ROL::Ptr< std::vector< Real > > &xvec, ROL::Vector< Real > &x) const
ROL::Ptr< BurgersFEM< Real > > fem_
H1VectorDual< Real > DualStateVector
H1VectorDual< Real > PrimalConstraintVector
void pruneLowerActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &x, Real eps=Real(0))
Set variables to zero if they correspond to the lower -active set.
Real apply(const ROL::Vector< Real > &x) const
Apply to a dual vector. This is equivalent to the call .
void scale(const Real alpha)
Compute where .
std::vector< Real >::size_type uint
void pruneLowerActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &g, const ROL::Vector< Real > &x, Real xeps=Real(0), Real geps=Real(0))
Set variables to zero if they correspond to the -binding set.
L2VectorPrimal(const ROL::Ptr< std::vector< Real > > &vec, const ROL::Ptr< BurgersFEM< Real > > &fem)
Real dot(const ROL::Vector< Real > &x) const override
Compute where .
Real dot(const ROL::Vector< Real > &x) const
Compute where .
H1VectorPrimal< Real > PrimalStateVector
void hessVec_22(ROL::Vector< Real > &hv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Constraint_BurgersControl(ROL::Ptr< BurgersFEM< Real > > &fem, bool useHessian=true)
Real value(const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Compute value.
void plus(const ROL::Vector< Real > &x)
Compute , where .
void gradient_1(ROL::Vector< Real > &g, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Compute gradient with respect to first component.
Real mesh_spacing(void) const
ROL::Ptr< ROL::Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
bool isFeasible(const ROL::Vector< Real > &x)
Check if the vector, v, is feasible.
int dimension() const
Return dimension of the vector space.
ROL::Ptr< std::vector< Real > > vec_
void test_inverse_mass(std::ostream &outStream=std::cout)
void apply_pde_jacobian(std::vector< Real > &jv, const std::vector< Real > &v, const std::vector< Real > &u, const std::vector< Real > &z) const
ROL::Ptr< BurgersFEM< Real > > fem_
const ROL::Ptr< const ROL::Vector< Real > > getLowerBound(void) const
Return the ref count pointer to the lower bound vector.
L2VectorBatchManager(const ROL::Ptr< const Teuchos::Comm< int > > &comm)
ROL::Ptr< ROL::Vector< Real > > basis(const int i) const
Return i-th basis vector.
void applyAdjointHessian_21(ROL::Vector< Real > &ahwv, const ROL::Vector< Real > &w, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Apply the simulation-space derivative of the adjoint of the constraint optimization-space Jacobian at...
void pruneUpperActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &x, Real eps=Real(0))
Set variables to zero if they correspond to the upper -active set.
ROL::Ptr< const std::vector< Real > > getVector() const
ROL::Ptr< ROL::Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
ROL::Ptr< const std::vector< Real > > getVector() const
ROL::Ptr< BurgersFEM< Real > > fem_
void set(const ROL::Vector< Real > &x)
Set where .
void scale(const Real alpha)
Compute where .
void compute_pde_jacobian(std::vector< Real > &dl, std::vector< Real > &d, std::vector< Real > &du, const std::vector< Real > &u) const
ROL::Ptr< ROL::Vector< Real > > basis(const int i) const
Return i-th basis vector.
H1VectorDual< Real > DualStateVector
void applyJacobian_1(ROL::Vector< Real > &jv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Apply the partial constraint Jacobian at , , to the vector .
Real compute_L2_dot(const std::vector< Real > &x, const std::vector< Real > &y) const
void hessVec_21(ROL::Vector< Real > &hv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
void linear_solve(std::vector< Real > &u, std::vector< Real > &dl, std::vector< Real > &d, std::vector< Real > &du, const std::vector< Real > &r, const bool transpose=false) const
std::vector< Real > x_lo_
ROL::Ptr< std::vector< Real > > getVector()
L2VectorPrimal< Real > PrimalControlVector
ROL::Ptr< ROL::Vector< Real > > u_
void apply_inverse_mass(std::vector< Real > &Mu, const std::vector< Real > &u) const
void pruneActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &g, const ROL::Vector< Real > &x, Real xeps=Real(0), Real geps=Real(0))
void axpy(std::vector< Real > &out, const Real a, const std::vector< Real > &x, const std::vector< Real > &y) const
void hessVec_12(ROL::Vector< Real > &hv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
ROL::Ptr< std::vector< Real > > getVector()
ROL::Ptr< BurgersFEM< Real > > fem_
ROL::Ptr< const std::vector< Real > > getVector() const
ROL::Ptr< BurgersFEM< Real > > fem_
void update(std::vector< Real > &u, const std::vector< Real > &s, const Real alpha=1.0) const
void pruneUpperActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &x, Real eps=Real(0))
Set variables to zero if they correspond to the upper -active set.
void scale(const Real alpha)
Compute where .
void projection(std::vector< Real > &x)
void gradient_2(ROL::Vector< Real > &g, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Compute gradient with respect to second component.
Real random(const ROL::Ptr< const Teuchos::Comm< int > > &comm)
void pruneActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &x, Real eps=Real(0))
void set(const ROL::Vector< Real > &x)
Set where .
ROL::Ptr< ROL::Vector< Real > > l_
Real apply(const ROL::Vector< Real > &x) const
Apply to a dual vector. This is equivalent to the call .
void apply_adjoint_control_hessian(std::vector< Real > &ahwv, const std::vector< Real > &w, const std::vector< Real > &v, const std::vector< Real > &u, const std::vector< Real > &z)
void applyInverseJacobian_1(ROL::Vector< Real > &ijv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Apply the inverse partial constraint Jacobian at , , to the vector .
Provides the interface to apply upper and lower bound constraints.
void compute_residual(std::vector< Real > &r, const std::vector< Real > &u, const std::vector< Real > &z) const
void projection(std::vector< Real > &x)
void cast_const_vector(ROL::Ptr< const std::vector< Real > > &xvec, const ROL::Vector< Real > &x) const
ROL::Ptr< ROL::Vector< Real > > l_
ROL::Ptr< ROL::Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
H1VectorBatchManager(const ROL::Ptr< const Teuchos::Comm< int > > &comm)
void set_problem_data(const Real nu, const Real f, const Real u0, const Real u1)
void hessVec_11(ROL::Vector< Real > &hv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Apply Hessian approximation to vector.
void applyJacobian_2(ROL::Vector< Real > &jv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Apply the partial constraint Jacobian at , , to the vector .
void apply_control_jacobian(std::vector< Real > &jv, const std::vector< Real > &v, const std::vector< Real > &u, const std::vector< Real > &z) const
Real dot(const ROL::Vector< Real > &x) const override
Compute where .
void test_inverse_H1(std::ostream &outStream=std::cout)
const ROL::Vector< Real > & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
void apply_adjoint_control_jacobian(std::vector< Real > &jv, const std::vector< Real > &v, const std::vector< Real > &u, const std::vector< Real > &z) const
Real dot(const ROL::Vector< Real > &x) const
Compute where .
void cast_vector(ROL::Ptr< std::vector< Real > > &xvec, ROL::Vector< Real > &x) const
L2VectorPrimal< Real > PrimalControlVector
void cast_const_vector(ROL::Ptr< const std::vector< Real > > &xvec, const ROL::Vector< Real > &x) const
H1VectorPrimal< Real > PrimalStateVector
std::vector< Real > x_lo_
void apply_inverse_pde_jacobian(std::vector< Real > &ijv, const std::vector< Real > &v, const std::vector< Real > &u, const std::vector< Real > &z) const
void axpy(std::vector< Real > &out, const Real a, const std::vector< Real > &x, const std::vector< Real > &y) const
L2VectorDual< Real > DualControlVector
void pruneLowerActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &g, const ROL::Vector< Real > &x, Real xeps=Real(0), Real geps=Real(0))
Set variables to zero if they correspond to the -binding set.
std::vector< int > indices_
ROL::Ptr< ROL::Vector< Real > > basis(const int i) const
Return i-th basis vector.
void applyAdjointHessian_22(ROL::Vector< Real > &ahwv, const ROL::Vector< Real > &w, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Apply the optimization-space derivative of the adjoint of the constraint optimization-space Jacobian ...
void apply_H1(std::vector< Real > &Mu, const std::vector< Real > &u) const
void plus(const ROL::Vector< Real > &x)
Compute , where .
const ROL::Ptr< const ROL::Vector< Real > > getLowerBound(void) const
Return the ref count pointer to the lower bound vector.
void sumAll(ROL::Vector< Real > &input, ROL::Vector< Real > &output)
void apply_inverse_H1(std::vector< Real > &Mu, const std::vector< Real > &u) const
void project(ROL::Vector< Real > &x)
Project optimization variables onto the bounds.
void apply_adjoint_pde_jacobian(std::vector< Real > &ajv, const std::vector< Real > &v, const std::vector< Real > &u, const std::vector< Real > &z) const
void set(const ROL::Vector< Real > &x)
Set where .
ROL::Ptr< BurgersFEM< Real > > fem_
ROL::Ptr< ROL::Vector< Real > > u_
Real apply(const ROL::Vector< Real > &x) const
Apply to a dual vector. This is equivalent to the call .
ROL::Ptr< ROL::Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
void apply_inverse_adjoint_pde_jacobian(std::vector< Real > &iajv, const std::vector< Real > &v, const std::vector< Real > &u, const std::vector< Real > &z) const
void apply_adjoint_control_pde_hessian(std::vector< Real > &ahwv, const std::vector< Real > &w, const std::vector< Real > &v, const std::vector< Real > &u, const std::vector< Real > &z)
void pruneUpperActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &g, const ROL::Vector< Real > &x, Real xeps=Real(0), Real geps=Real(0))
Set variables to zero if they correspond to the upper -binding set.
Objective_BurgersControl(const ROL::Ptr< BurgersFEM< Real > > &fem, Real x=0.0)
void pruneActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &x, Real eps=Real(0))
void applyAdjointJacobian_2(ROL::Vector< Real > &jv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Apply the adjoint of the partial constraint Jacobian at , , to vector . This is the primary interface...
int dimension() const
Return dimension of the vector space.
Defines the constraint operator interface for simulation-based optimization.
L2VectorDual< Real > DualControlVector
ROL::Ptr< H1VectorPrimal< Real > > dual_vec_
const ROL::Vector< Real > & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
void value(ROL::Vector< Real > &c, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Evaluate the constraint operator at .
void sumAll(ROL::Vector< Real > &input, ROL::Vector< Real > &output)
void set(const ROL::Vector< Real > &x)
Set where .
Real norm() const
Returns where .
void scale(std::vector< Real > &u, const Real alpha=0.0) const
const ROL::Vector< Real > & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
ROL::Ptr< H1VectorDual< Real > > dual_vec_
const ROL::Ptr< const ROL::Vector< Real > > getUpperBound(void) const
Return the ref count pointer to the upper bound vector.
Real dot(const ROL::Vector< Real > &x) const override
Compute where .
ROL::Ptr< std::vector< Real > > getVector()
std::vector< Real >::size_type uint
H1BoundConstraint(std::vector< Real > &l, std::vector< Real > &u, const ROL::Ptr< BurgersFEM< Real > > &fem, Real scale=1.0)
void scale(const Real alpha)
Compute where .
int dimension() const
Return dimension of the vector space.
const ROL::Vector< Real > & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
Real compute_L2_norm(const std::vector< Real > &r) const
ROL::Ptr< L2VectorPrimal< Real > > dual_vec_