20 #include "ROL_TeuchosBatchManager.hpp"
22 #include "Teuchos_LAPACK.hpp"
50 void update(std::vector<Real> &u,
const std::vector<Real> &s,
const Real alpha=1.0)
const {
51 for (
unsigned i=0; i<u.size(); i++) {
56 void axpy(std::vector<Real> &out,
const Real a,
const std::vector<Real> &x,
const std::vector<Real> &y)
const {
57 for (
unsigned i=0; i < x.size(); i++) {
58 out[i] = a*x[i] + y[i];
62 void scale(std::vector<Real> &u,
const Real alpha=0.0)
const {
63 for (
unsigned i=0; i<u.size(); i++) {
68 void linear_solve(std::vector<Real> &u, std::vector<Real> &dl, std::vector<Real> &d, std::vector<Real> &du,
69 const std::vector<Real> &r,
const bool transpose =
false)
const {
70 if ( r.size() == 1 ) {
71 u.resize(1,r[0]/d[0]);
73 else if ( r.size() == 2 ) {
75 Real det = d[0]*d[1] - dl[0]*du[0];
76 u[0] = (d[1]*r[0] - du[0]*r[1])/det;
77 u[1] = (d[0]*r[1] - dl[0]*r[0])/det;
80 u.assign(r.begin(),r.end());
82 Teuchos::LAPACK<int,Real> lp;
83 std::vector<Real> du2(r.size()-2,0.0);
84 std::vector<int> ipiv(r.size(),0);
89 lp.GTTRF(dim,&dl[0],&d[0],&du[0],&du2[0],&ipiv[0],&info);
94 lp.GTTRS(trans,dim,nhrs,&dl[0],&d[0],&du[0],&du2[0],&ipiv[0],&u[0],ldb,&info);
99 BurgersFEM(
int nx = 128, Real nl = 1.0, Real cH1 = 1.0, Real cL2 = 1.0)
103 nu_ = std::pow(10.0,nu-2.0);
121 Real
compute_L2_dot(
const std::vector<Real> &x,
const std::vector<Real> &y)
const {
123 Real c = (((int)x.size()==
nx_) ? 4.0 : 2.0);
124 for (
unsigned i=0; i<x.size(); i++) {
126 ip +=
dx_/6.0*(c*x[i] + x[i+1])*y[i];
128 else if ( i == x.size()-1 ) {
129 ip +=
dx_/6.0*(x[i-1] + c*x[i])*y[i];
132 ip +=
dx_/6.0*(x[i-1] + 4.0*x[i] + x[i+1])*y[i];
144 void apply_mass(std::vector<Real> &Mu,
const std::vector<Real> &u )
const {
145 Mu.resize(u.size(),0.0);
146 Real c = (((int)u.size()==
nx_) ? 4.0 : 2.0);
147 for (
unsigned i=0; i<u.size(); i++) {
149 Mu[i] =
dx_/6.0*(c*u[i] + u[i+1]);
151 else if ( i == u.size()-1 ) {
152 Mu[i] =
dx_/6.0*(u[i-1] + c*u[i]);
155 Mu[i] =
dx_/6.0*(u[i-1] + 4.0*u[i] + u[i+1]);
162 unsigned nx = u.size();
164 std::vector<Real> dl(nx-1,
dx_/6.0);
165 std::vector<Real> d(nx,2.0*
dx_/3.0);
166 std::vector<Real> du(nx-1,
dx_/6.0);
167 if ( (
int)nx !=
nx_ ) {
176 std::vector<Real> u(
nx_,0.0), Mu(
nx_,0.0), iMMu(
nx_,0.0), diff(
nx_,0.0);
177 for (
int i = 0; i <
nx_; i++) {
178 u[i] = 2.0*(Real)rand()/(Real)RAND_MAX - 1.0;
182 axpy(diff,-1.0,iMMu,u);
185 outStream <<
"Test Inverse State Mass Matrix\n";
186 outStream <<
" ||u - inv(M)Mu|| = " << error <<
"\n";
187 outStream <<
" ||u|| = " << normu <<
"\n";
188 outStream <<
" Relative Error = " << error/normu <<
"\n";
191 u.resize(nx_+2,0.0); Mu.resize(nx_+2,0.0); iMMu.resize(nx_+2,0.0); diff.resize(nx_+2,0.0);
192 for (
int i = 0; i < nx_+2; i++) {
193 u[i] = 2.0*(Real)rand()/(Real)RAND_MAX - 1.0;
197 axpy(diff,-1.0,iMMu,u);
200 outStream <<
"Test Inverse Control Mass Matrix\n";
201 outStream <<
" ||z - inv(M)Mz|| = " << error <<
"\n";
202 outStream <<
" ||z|| = " << normu <<
"\n";
203 outStream <<
" Relative Error = " << error/normu <<
"\n";
211 Real
compute_H1_dot(
const std::vector<Real> &x,
const std::vector<Real> &y)
const {
213 for (
int i=0; i<
nx_; i++) {
215 ip +=
cL2_*
dx_/6.0*(4.0*x[i] + x[i+1])*y[i];
216 ip +=
cH1_*(2.0*x[i] - x[i+1])/
dx_*y[i];
218 else if ( i == nx_-1 ) {
219 ip +=
cL2_*
dx_/6.0*(x[i-1] + 4.0*x[i])*y[i];
220 ip +=
cH1_*(2.0*x[i] - x[i-1])/
dx_*y[i];
223 ip +=
cL2_*
dx_/6.0*(x[i-1] + 4.0*x[i] + x[i+1])*y[i];
224 ip +=
cH1_*(2.0*x[i] - x[i-1] - x[i+1])/
dx_*y[i];
236 void apply_H1(std::vector<Real> &Mu,
const std::vector<Real> &u )
const {
238 for (
int i=0; i<
nx_; i++) {
240 Mu[i] =
cL2_*
dx_/6.0*(4.0*u[i] + u[i+1])
241 +
cH1_*(2.0*u[i] - u[i+1])/
dx_;
243 else if ( i == nx_-1 ) {
244 Mu[i] =
cL2_*
dx_/6.0*(u[i-1] + 4.0*u[i])
245 +
cH1_*(2.0*u[i] - u[i-1])/
dx_;
248 Mu[i] =
cL2_*
dx_/6.0*(u[i-1] + 4.0*u[i] + u[i+1])
249 +
cH1_*(2.0*u[i] - u[i-1] - u[i+1])/
dx_;
264 std::vector<Real> u(
nx_,0.0), Mu(
nx_,0.0), iMMu(
nx_,0.0), diff(
nx_,0.0);
265 for (
int i = 0; i <
nx_; i++) {
266 u[i] = 2.0*(Real)rand()/(Real)RAND_MAX - 1.0;
270 axpy(diff,-1.0,iMMu,u);
273 outStream <<
"Test Inverse State H1 Matrix\n";
274 outStream <<
" ||u - inv(M)Mu|| = " << error <<
"\n";
275 outStream <<
" ||u|| = " << normu <<
"\n";
276 outStream <<
" Relative Error = " << error/normu <<
"\n";
285 const std::vector<Real> &z)
const {
288 for (
int i=0; i<
nx_; i++) {
291 r[i] =
nu_/
dx_*(2.0*u[i]-u[i+1]);
294 r[i] =
nu_/
dx_*(2.0*u[i]-u[i-1]);
297 r[i] =
nu_/
dx_*(2.0*u[i]-u[i-1]-u[i+1]);
301 r[i] +=
nl_*u[i+1]*(u[i]+u[i+1])/6.0;
304 r[i] -=
nl_*u[i-1]*(u[i-1]+u[i])/6.0;
307 r[i] -=
dx_/6.0*(z[i]+4.0*z[i+1]+z[i+2]);
321 std::vector<Real> &d,
322 std::vector<Real> &du,
323 const std::vector<Real> &u)
const {
332 for (
int i=0; i<
nx_; i++) {
334 dl[i] +=
nl_*(-2.0*u[i]-u[i+1])/6.0;
335 d[i] +=
nl_*u[i+1]/6.0;
338 d[i] -=
nl_*u[i-1]/6.0;
339 du[i-1] +=
nl_*(u[i-1]+2.0*u[i])/6.0;
349 const std::vector<Real> &v,
350 const std::vector<Real> &u,
351 const std::vector<Real> &z)
const {
353 for (
int i = 0; i <
nx_; i++) {
356 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]);
359 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]);
362 jv[ 0] -=
nl_*
u0_/6.0*v[0];
363 jv[nx_-1] +=
nl_*
u1_/6.0*v[nx_-1];
368 const std::vector<Real> &v,
369 const std::vector<Real> &u,
370 const std::vector<Real> &z)
const {
372 std::vector<Real> d(
nx_,0.0);
373 std::vector<Real> dl(
nx_-1,0.0);
374 std::vector<Real> du(
nx_-1,0.0);
382 const std::vector<Real> &v,
383 const std::vector<Real> &u,
384 const std::vector<Real> &z)
const {
386 for (
int i = 0; i <
nx_; i++) {
387 ajv[i] =
nu_/
dx_*2.0*v[i];
389 ajv[i] += -
nu_/
dx_*v[i-1]-
nl_*(u[i-1]/6.0*v[i]
390 -(u[i-1]+2.0*u[i])/6.0*v[i-1]);
393 ajv[i] += -
nu_/
dx_*v[i+1]+
nl_*(u[i+1]/6.0*v[i]
394 -(u[i+1]+2.0*u[i])/6.0*v[i+1]);
397 ajv[ 0] -=
nl_*
u0_/6.0*v[0];
398 ajv[nx_-1] +=
nl_*
u1_/6.0*v[nx_-1];
403 const std::vector<Real> &v,
404 const std::vector<Real> &u,
405 const std::vector<Real> &z)
const {
407 std::vector<Real> d(
nx_,0.0);
408 std::vector<Real> du(
nx_-1,0.0);
409 std::vector<Real> dl(
nx_-1,0.0);
420 const std::vector<Real> &v,
421 const std::vector<Real> &u,
422 const std::vector<Real> &z)
const {
423 for (
int i=0; i<
nx_; i++) {
425 jv[i] = -
dx_/6.0*(v[i]+4.0*v[i+1]+v[i+2]);
431 const std::vector<Real> &v,
432 const std::vector<Real> &u,
433 const std::vector<Real> &z)
const {
434 for (
int i=0; i<
nx_+2; i++) {
436 jv[i] = -
dx_/6.0*v[i];
439 jv[i] = -
dx_/6.0*(4.0*v[i-1]+v[i]);
441 else if ( i == nx_ ) {
442 jv[i] = -
dx_/6.0*(4.0*v[i-1]+v[i-2]);
444 else if ( i == nx_+1 ) {
445 jv[i] = -
dx_/6.0*v[i-2];
448 jv[i] = -
dx_/6.0*(v[i-2]+4.0*v[i-1]+v[i]);
457 const std::vector<Real> &w,
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++) {
465 ahwv[i] += (w[i]*v[i+1] - w[i+1]*(2.0*v[i]+v[i+1]))/6.0;
468 ahwv[i] += (w[i-1]*(v[i-1]+2.0*v[i]) - w[i]*v[i-1])/6.0;
474 const std::vector<Real> &w,
475 const std::vector<Real> &v,
476 const std::vector<Real> &u,
477 const std::vector<Real> &z) {
478 ahwv.assign(u.size(),0.0);
481 const std::vector<Real> &w,
482 const std::vector<Real> &v,
483 const std::vector<Real> &u,
484 const std::vector<Real> &z) {
485 ahwv.assign(z.size(),0.0);
488 const std::vector<Real> &w,
489 const std::vector<Real> &v,
490 const std::vector<Real> &u,
491 const std::vector<Real> &z) {
492 ahwv.assign(z.size(),0.0);
499 ROL::Ptr<BurgersFEM<Real> >
fem_;
500 mutable ROL::Ptr<L2VectorDual<Real> >
dual_vec_;
511 const std::vector<Real>& xval = *ex.
getVector();
513 return fem_->compute_L2_dot(xval,yval);
516 ROL::Ptr<ROL::Vector<Real> >
clone()
const {
517 return ROL::makePtr<L2VectorPrimal>(
523 dual_vec_ = ROL::makePtr<L2VectorDual<Real>>(
536 ROL::Ptr<BurgersFEM<Real> >
fem_;
548 const std::vector<Real>& xval = *ex.
getVector();
551 std::vector<Real> Mx(dimension,0.0);
552 fem_->apply_inverse_mass(Mx,xval);
554 for (
unsigned i = 0; i <
dimension; i++) {
555 val += Mx[i]*yval[i];
560 ROL::Ptr<ROL::Vector<Real> >
clone()
const {
561 return ROL::makePtr<L2VectorDual>(
567 prim_vec_ = ROL::makePtr<L2VectorPrimal<Real>>(
580 ROL::Ptr<BurgersFEM<Real> >
fem_;
581 mutable ROL::Ptr<H1VectorDual<Real> >
dual_vec_;
592 const std::vector<Real>& xval = *ex.
getVector();
594 return fem_->compute_H1_dot(xval,yval);
597 ROL::Ptr<ROL::Vector<Real> >
clone()
const {
598 return ROL::makePtr<H1VectorPrimal>(
604 dual_vec_ = ROL::makePtr<H1VectorDual<Real>>(
617 ROL::Ptr<BurgersFEM<Real> >
fem_;
629 const std::vector<Real>& xval = *ex.
getVector();
632 std::vector<Real> Mx(dimension,0.0);
633 fem_->apply_inverse_H1(Mx,xval);
635 for (
unsigned i = 0; i <
dimension; i++) {
636 val += Mx[i]*yval[i];
641 ROL::Ptr<ROL::Vector<Real> >
clone()
const {
642 return ROL::makePtr<H1VectorDual>(
648 prim_vec_ = ROL::makePtr<H1VectorPrimal<Real>>(
671 ROL::Ptr<BurgersFEM<Real> >
fem_;
680 ROL::Ptr<std::vector<Real> > cp =
682 ROL::Ptr<const std::vector<Real> > up =
684 ROL::Ptr<const std::vector<Real> > zp =
687 const std::vector<Real> param
689 fem_->set_problem_data(param[0],param[1],param[2],param[3]);
691 fem_->compute_residual(*cp,*up,*zp);
696 ROL::Ptr<std::vector<Real> > jvp =
698 ROL::Ptr<const std::vector<Real> > vp =
700 ROL::Ptr<const std::vector<Real> > up =
702 ROL::Ptr<const std::vector<Real> > zp =
705 const std::vector<Real> param
707 fem_->set_problem_data(param[0],param[1],param[2],param[3]);
709 fem_->apply_pde_jacobian(*jvp,*vp,*up,*zp);
714 ROL::Ptr<std::vector<Real> > jvp =
716 ROL::Ptr<const std::vector<Real> > vp =
718 ROL::Ptr<const std::vector<Real> > up =
720 ROL::Ptr<const std::vector<Real> > zp =
723 const std::vector<Real> param
725 fem_->set_problem_data(param[0],param[1],param[2],param[3]);
727 fem_->apply_control_jacobian(*jvp,*vp,*up,*zp);
732 ROL::Ptr<std::vector<Real> > ijvp =
734 ROL::Ptr<const std::vector<Real> > vp =
736 ROL::Ptr<const std::vector<Real> > up =
738 ROL::Ptr<const std::vector<Real> > zp =
741 const std::vector<Real> param
743 fem_->set_problem_data(param[0],param[1],param[2],param[3]);
745 fem_->apply_inverse_pde_jacobian(*ijvp,*vp,*up,*zp);
750 ROL::Ptr<std::vector<Real> > jvp =
752 ROL::Ptr<const std::vector<Real> > vp =
754 ROL::Ptr<const std::vector<Real> > up =
756 ROL::Ptr<const std::vector<Real> > zp =
759 const std::vector<Real> param
761 fem_->set_problem_data(param[0],param[1],param[2],param[3]);
763 fem_->apply_adjoint_pde_jacobian(*jvp,*vp,*up,*zp);
768 ROL::Ptr<std::vector<Real> > jvp =
770 ROL::Ptr<const std::vector<Real> > vp =
772 ROL::Ptr<const std::vector<Real> > up =
774 ROL::Ptr<const std::vector<Real> > zp =
777 const std::vector<Real> param
779 fem_->set_problem_data(param[0],param[1],param[2],param[3]);
781 fem_->apply_adjoint_control_jacobian(*jvp,*vp,*up,*zp);
786 ROL::Ptr<std::vector<Real> > iajvp =
788 ROL::Ptr<const std::vector<Real> > vp =
790 ROL::Ptr<const std::vector<Real> > up =
792 ROL::Ptr<const std::vector<Real> > zp =
795 const std::vector<Real> param
797 fem_->set_problem_data(param[0],param[1],param[2],param[3]);
799 fem_->apply_inverse_adjoint_pde_jacobian(*iajvp,*vp,*up,*zp);
805 ROL::Ptr<std::vector<Real> > ahwvp =
807 ROL::Ptr<const std::vector<Real> > wp =
809 ROL::Ptr<const std::vector<Real> > vp =
811 ROL::Ptr<const std::vector<Real> > up =
813 ROL::Ptr<const std::vector<Real> > zp =
816 const std::vector<Real> param
818 fem_->set_problem_data(param[0],param[1],param[2],param[3]);
820 fem_->apply_adjoint_pde_hessian(*ahwvp,*wp,*vp,*up,*zp);
830 ROL::Ptr<std::vector<Real> > ahwvp =
832 ROL::Ptr<const std::vector<Real> > wp =
834 ROL::Ptr<const std::vector<Real> > vp =
836 ROL::Ptr<const std::vector<Real> > up =
838 ROL::Ptr<const std::vector<Real> > zp =
841 const std::vector<Real> param
843 fem_->set_problem_data(param[0],param[1],param[2],param[3]);
845 fem_->apply_adjoint_control_pde_hessian(*ahwvp,*wp,*vp,*up,*zp);
854 ROL::Ptr<std::vector<Real> > ahwvp =
856 ROL::Ptr<const std::vector<Real> > wp =
858 ROL::Ptr<const std::vector<Real> > vp =
860 ROL::Ptr<const std::vector<Real> > up =
862 ROL::Ptr<const std::vector<Real> > zp =
865 const std::vector<Real> param
867 fem_->set_problem_data(param[0],param[1],param[2],param[3]);
869 fem_->apply_adjoint_pde_control_hessian(*ahwvp,*wp,*vp,*up,*zp);
878 ROL::Ptr<std::vector<Real> > ahwvp =
880 ROL::Ptr<const std::vector<Real> > wp =
882 ROL::Ptr<const std::vector<Real> > vp =
884 ROL::Ptr<const std::vector<Real> > up =
886 ROL::Ptr<const std::vector<Real> > zp =
889 const std::vector<Real> param
891 fem_->set_problem_data(param[0],param[1],param[2],param[3]);
893 fem_->apply_adjoint_control_hessian(*ahwvp,*wp,*vp,*up,*zp);
912 ROL::Ptr<BurgersFEM<Real> >
fem_;
913 ROL::Ptr<ROL::Vector<Real> >
ud_;
914 ROL::Ptr<ROL::Vector<Real> >
diff_;
926 ROL::Ptr<const std::vector<Real> > up =
928 ROL::Ptr<const std::vector<Real> > zp =
930 ROL::Ptr<const std::vector<Real> > udp =
933 std::vector<Real> diff(udp->size(),0.0);
934 for (
unsigned i = 0; i < udp->size(); i++) {
935 diff[i] = (*up)[i] - (*udp)[i];
937 return 0.5*(
fem_->compute_L2_dot(diff,diff) +
alpha_*
fem_->compute_L2_dot(*zp,*zp));
941 ROL::Ptr<std::vector<Real> > gp =
943 ROL::Ptr<const std::vector<Real> > up =
945 ROL::Ptr<const std::vector<Real> > udp =
948 std::vector<Real> diff(udp->size(),0.0);
949 for (
unsigned i = 0; i < udp->size(); i++) {
950 diff[i] = (*up)[i] - (*udp)[i];
952 fem_->apply_mass(*gp,diff);
956 ROL::Ptr<std::vector<Real> > gp =
958 ROL::Ptr<const std::vector<Real> > zp =
961 fem_->apply_mass(*gp,*zp);
962 for (
unsigned i = 0; i < zp->size(); i++) {
969 ROL::Ptr<std::vector<Real> > hvp =
971 ROL::Ptr<const std::vector<Real> > vp =
974 fem_->apply_mass(*hvp,*vp);
989 ROL::Ptr<std::vector<Real> > hvp =
991 ROL::Ptr<const std::vector<Real> > vp =
994 fem_->apply_mass(*hvp,*vp);
995 for (
unsigned i = 0; i < vp->size(); i++) {
1001 template<
class Real>
1005 std::vector<Real>
x_lo_;
1006 std::vector<Real>
x_up_;
1009 ROL::Ptr<BurgersFEM<Real> >
fem_;
1010 ROL::Ptr<ROL::Vector<Real> >
l_;
1011 ROL::Ptr<ROL::Vector<Real> >
u_;
1018 catch (std::exception &e) {
1028 catch (std::exception &e) {
1033 void axpy(std::vector<Real> &out,
const Real a,
1034 const std::vector<Real> &x,
const std::vector<Real> &y)
const{
1035 out.resize(
dim_,0.0);
1036 for (
unsigned i = 0; i <
dim_; i++) {
1037 out[i] = a*x[i] + y[i];
1042 for (
int i = 0; i <
dim_; i++ ) {
1043 x[i] = std::max(
x_lo_[i],std::min(
x_up_[i],x[i]));
1052 for (
int i = 0; i <
dim_; i++ ) {
1061 l_ = ROL::makePtr<L2VectorPrimal<Real>>(
1062 ROL::makePtr<std::vector<Real>>(l), fem);
1063 u_ = ROL::makePtr<L2VectorPrimal<Real>>(
1064 ROL::makePtr<std::vector<Real>>(u), fem);
1071 for (
int i = 0; i <
dim_; i++ ) {
1072 if ( (*ex)[i] >=
x_lo_[i] && (*ex)[i] <=
x_up_[i] ) { cnt *= 1; }
1075 if ( cnt == 0 ) { val =
false; }
1080 ROL::Ptr<std::vector<Real> > ex;
cast_vector(ex,x);
1086 ROL::Ptr<std::vector<Real> > ev;
cast_vector(ev,v);
1088 for (
int i = 0; i <
dim_; i++ ) {
1089 if ( ((*ex)[i] <=
x_lo_[i]+epsn) ) {
1097 ROL::Ptr<std::vector<Real> > ev;
cast_vector(ev,v);
1099 for (
int i = 0; i <
dim_; i++ ) {
1100 if ( ((*ex)[i] >=
x_up_[i]-epsn) ) {
1108 ROL::Ptr<std::vector<Real> > ev;
cast_vector(ev,v);
1110 for (
int i = 0; i <
dim_; i++ ) {
1111 if ( ((*ex)[i] <=
x_lo_[i]+epsn) ||
1112 ((*ex)[i] >=
x_up_[i]-epsn) ) {
1121 ROL::Ptr<std::vector<Real> > ev;
cast_vector(ev,v);
1123 for (
int i = 0; i <
dim_; i++ ) {
1124 if ( ((*ex)[i] <=
x_lo_[i]+epsn && (*eg)[i] > geps) ) {
1133 ROL::Ptr<std::vector<Real> > ev;
cast_vector(ev,v);
1135 for (
int i = 0; i <
dim_; i++ ) {
1136 if ( ((*ex)[i] >=
x_up_[i]-epsn && (*eg)[i] < -geps) ) {
1145 ROL::Ptr<std::vector<Real> > ev;
cast_vector(ev,v);
1147 for (
int i = 0; i <
dim_; i++ ) {
1148 if ( ((*ex)[i] <=
x_lo_[i]+epsn && (*eg)[i] > geps) ||
1149 ((*ex)[i] >=
x_up_[i]-epsn && (*eg)[i] < -geps) ) {
1164 template<
class Real>
1168 std::vector<Real>
x_lo_;
1169 std::vector<Real>
x_up_;
1172 ROL::Ptr<BurgersFEM<Real> >
fem_;
1173 ROL::Ptr<ROL::Vector<Real> >
l_;
1174 ROL::Ptr<ROL::Vector<Real> >
u_;
1179 xvec = ROL::constPtrCast<std::vector<Real> >(
1182 catch (std::exception &e) {
1183 xvec = ROL::constPtrCast<std::vector<Real> >(
1193 catch (std::exception &e) {
1198 void axpy(std::vector<Real> &out,
const Real a,
1199 const std::vector<Real> &x,
const std::vector<Real> &y)
const{
1200 out.resize(
dim_,0.0);
1201 for (
unsigned i = 0; i <
dim_; i++) {
1202 out[i] = a*x[i] + y[i];
1207 for (
int i = 0; i <
dim_; i++ ) {
1208 x[i] = std::max(
x_lo_[i],std::min(
x_up_[i],x[i]));
1217 for (
int i = 0; i <
dim_; i++ ) {
1226 l_ = ROL::makePtr<H1VectorPrimal<Real>>(
1227 ROL::makePtr<std::vector<Real>>(l), fem);
1228 u_ = ROL::makePtr<H1VectorPrimal<Real>>(
1229 ROL::makePtr<std::vector<Real>>(u), fem);
1236 for (
int i = 0; i <
dim_; i++ ) {
1237 if ( (*ex)[i] >=
x_lo_[i] && (*ex)[i] <=
x_up_[i] ) { cnt *= 1; }
1240 if ( cnt == 0 ) { val =
false; }
1245 ROL::Ptr<std::vector<Real> > ex;
cast_vector(ex,x);
1251 ROL::Ptr<std::vector<Real> > ev;
cast_vector(ev,v);
1253 for (
int i = 0; i <
dim_; i++ ) {
1254 if ( ((*ex)[i] <=
x_lo_[i]+epsn) ) {
1262 ROL::Ptr<std::vector<Real> > ev;
cast_vector(ev,v);
1264 for (
int i = 0; i <
dim_; i++ ) {
1265 if ( ((*ex)[i] >=
x_up_[i]-epsn) ) {
1273 ROL::Ptr<std::vector<Real> > ev;
cast_vector(ev,v);
1275 for (
int i = 0; i <
dim_; i++ ) {
1276 if ( ((*ex)[i] <=
x_lo_[i]+epsn) ||
1277 ((*ex)[i] >=
x_up_[i]-epsn) ) {
1286 ROL::Ptr<std::vector<Real> > ev;
cast_vector(ev,v);
1288 for (
int i = 0; i <
dim_; i++ ) {
1289 if ( ((*ex)[i] <=
x_lo_[i]+epsn && (*eg)[i] > geps) ) {
1298 ROL::Ptr<std::vector<Real> > ev;
cast_vector(ev,v);
1300 for (
int i = 0; i <
dim_; i++ ) {
1301 if ( ((*ex)[i] >=
x_up_[i]-epsn && (*eg)[i] < -geps) ) {
1310 ROL::Ptr<std::vector<Real> > ev;
cast_vector(ev,v);
1312 for (
int i = 0; i <
dim_; i++ ) {
1313 if ( ((*ex)[i] <=
x_lo_[i]+epsn && (*eg)[i] > geps) ||
1314 ((*ex)[i] >=
x_up_[i]-epsn && (*eg)[i] < -geps) ) {
1329 template<
class Real,
class Ordinal>
1337 catch (std::exception &e) {
1344 : ROL::TeuchosBatchManager<Real,Ordinal>(comm) {}
1346 ROL::Ptr<std::vector<Real> > input_ptr;
1348 int dim_i = input_ptr->size();
1349 ROL::Ptr<std::vector<Real> > output_ptr;
1351 int dim_o = output_ptr->size();
1352 if ( dim_i != dim_o ) {
1353 std::cout <<
"L2VectorBatchManager: DIMENSION MISMATCH ON RANK "
1354 << ROL::TeuchosBatchManager<Real,Ordinal>::batchID() <<
"\n";
1357 ROL::TeuchosBatchManager<Real,Ordinal>::sumAll(&(*input_ptr)[0],&(*output_ptr)[0],dim_i);
1362 template<
class Real,
class Ordinal>
1370 catch (std::exception &e) {
1377 : ROL::TeuchosBatchManager<Real,Ordinal>(comm) {}
1379 ROL::Ptr<std::vector<Real> > input_ptr;
1381 int dim_i = input_ptr->size();
1382 ROL::Ptr<std::vector<Real> > output_ptr;
1384 int dim_o = output_ptr->size();
1385 if ( dim_i != dim_o ) {
1386 std::cout <<
"H1VectorBatchManager: DIMENSION MISMATCH ON RANK "
1387 << ROL::TeuchosBatchManager<Real,Ordinal>::batchID() <<
"\n";
1390 ROL::TeuchosBatchManager<Real,Ordinal>::sumAll(&(*input_ptr)[0],&(*output_ptr)[0],dim_i);
1395 template<
class Real>
1396 Real
random(
const ROL::Ptr<
const Teuchos::Comm<int> > &comm) {
1398 if ( Teuchos::rank<int>(*comm)==0 ) {
1399 val = (Real)rand()/(Real)RAND_MAX;
1401 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)
Provides the interface to evaluate simulation-based objective functions.
std::vector< Real > x_up_
ROL::Ptr< ROL::Vector< Real > > diff_
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...
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
void cast_vector(ROL::Ptr< std::vector< Real > > &xvec, ROL::Vector< Real > &x) const
ROL::Ptr< L2VectorPrimal< Real > > prim_vec_
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< H1VectorPrimal< Real > > prim_vec_
ROL::Ptr< L2VectorDual< Real > > dual_vec_
StdVector(const Ptr< std::vector< Real >> &std_vec)
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...
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)
Contains definitions of custom data types in ROL.
Ptr< const std::vector< Element > > getVector() 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
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.
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)
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 ...
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.
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.
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
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 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.
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.
ROL::Ptr< ROL::Vector< Real > > ud_
L2VectorBatchManager(const ROL::Ptr< const Teuchos::Comm< int > > &comm)
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...
Provides the ROL::Vector interface for scalar values, to be used, for example, with scalar constraint...
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 compute_pde_jacobian(std::vector< Real > &dl, std::vector< Real > &d, std::vector< Real > &du, const std::vector< Real > &u) const
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_
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< 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 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))
int dimension() const
Return dimension of the vector space.
ROL::Ptr< ROL::Vector< Real > > l_
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
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.
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
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
ROL::Ptr< BurgersFEM< Real > > fem_
ROL::Ptr< ROL::Vector< Real > > u_
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.
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
const ROL::Vector< Real > & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
bool isPrimalInitialized_
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)
bool isPrimalInitialized_
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_
Objective_BurgersControl(const ROL::Ptr< BurgersFEM< Real > > &fem, const ROL::Ptr< ROL::Vector< Real > > &ud, Real alpha=1.e-4)
const ROL::Ptr< const ROL::Vector< Real > > getUpperBound(void) const
Return the ref count pointer to the upper bound vector.
H1BoundConstraint(std::vector< Real > &l, std::vector< Real > &u, const ROL::Ptr< BurgersFEM< Real > > &fem, Real scale=1.0)
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