51 #include "Teuchos_oblackholestream.hpp" 52 #include "Teuchos_GlobalMPISession.hpp" 53 #include "Teuchos_XMLParameterListHelpers.hpp" 54 #include "Teuchos_LAPACK.hpp" 78 case 1: val = ((x<0.5) ? 1.0 : 0.0);
break;
79 case 2: val = 1.0;
break;
80 case 3: val = std::abs(std::sin(8.0*M_PI*x));
break;
81 case 4: val = std::exp(-0.5*(x-0.5)*(x-0.5));
break;
87 return std::sqrt(this->
dot(r,r));
90 Real
dot(
const std::vector<Real> &x,
const std::vector<Real> &y) {
92 Real c = (((int)x.size()==this->
nx_) ? 4.0 : 2.0);
93 for (
unsigned i=0; i<x.size(); i++) {
95 ip += this->dx_/6.0*(c*x[i] + x[i+1])*y[i];
97 else if ( i == x.size()-1 ) {
98 ip += this->dx_/6.0*(x[i-1] + c*x[i])*y[i];
101 ip += this->dx_/6.0*(x[i-1] + 4.0*x[i] + x[i+1])*y[i];
107 void update(std::vector<Real> &u,
const std::vector<Real> &s,
const Real alpha=1.0) {
108 for (
unsigned i=0; i<u.size(); i++) {
113 void scale(std::vector<Real> &u,
const Real alpha=0.0) {
114 for (
unsigned i=0; i<u.size(); i++) {
119 void apply_mass(std::vector<Real> &Mu,
const std::vector<Real> &u ) {
120 Mu.resize(u.size(),0.0);
121 Real c = (((int)u.size()==this->
nx_) ? 4.0 : 2.0);
122 for (
unsigned i=0; i<u.size(); i++) {
124 Mu[i] = this->dx_/6.0*(c*u[i] + u[i+1]);
126 else if ( i == u.size()-1 ) {
127 Mu[i] = this->dx_/6.0*(u[i-1] + c*u[i]);
130 Mu[i] = this->dx_/6.0*(u[i-1] + 4.0*u[i] + u[i+1]);
136 const std::vector<Real> &z,
const std::vector<Real> ¶m) {
138 r.resize(this->nx_,0.0);
139 Real nu = std::pow(10.0,param[0]-2.0);
140 Real f = param[1]/100.0;
141 Real u0 = 1.0+param[2]/1000.0;
142 Real u1 = param[3]/1000.0;
143 for (
int i=0; i<this->
nx_; i++) {
146 r[i] = nu/this->dx_*(2.0*u[i]-u[i+1]);
148 else if (i==this->nx_-1) {
149 r[i] = nu/this->dx_*(2.0*u[i]-u[i-1]);
152 r[i] = nu/this->dx_*(2.0*u[i]-u[i-1]-u[i+1]);
156 r[i] += u[i+1]*(u[i]+u[i+1])/6.0;
159 r[i] -= u[i-1]*(u[i-1]+u[i])/6.0;
162 r[i] -= this->dx_/6.0*(z[i]+4.0*z[i+1]+z[i+2]);
167 r[0] -= u0*u[ 0]/6.0 + u0*u0/6.0 + nu*u0/this->
dx_;
168 r[this->nx_-1] += u1*u[this->nx_-1]/6.0 + u1*u1/6.0 - nu*u1/this->
dx_;
172 const std::vector<Real> &u,
const std::vector<Real> ¶m) {
173 Real nu = std::pow(10.0,param[0]-2.0);
174 Real u0 = 1.0+param[2]/1000.0;
175 Real u1 = param[3]/1000.0;
178 d.resize(this->nx_,nu*2.0/this->dx_);
180 dl.resize(this->nx_-1,-nu/this->dx_);
182 du.resize(this->nx_-1,-nu/this->dx_);
184 for (
int i=0; i<this->
nx_; i++) {
186 dl[i] += (-2.0*u[i]-u[i+1])/6.0;
191 du[i-1] += (u[i-1]+2.0*u[i])/6.0;
196 d[this->nx_-1] += u1/6.0;
199 void add_pde_hessian(std::vector<Real> &r,
const std::vector<Real> &u,
const std::vector<Real> &p,
200 const std::vector<Real> &s, Real alpha = 1.0) {
201 for (
int i=0; i<this->
nx_; i++) {
204 r[i] += alpha*(p[i]*s[i+1] - p[i+1]*(2.0*s[i]+s[i+1]))/6.0;
207 r[i] += alpha*(p[i-1]*(s[i-1]+2.0*s[i]) - p[i]*s[i-1])/6.0;
212 void linear_solve(std::vector<Real> &u, std::vector<Real> &dl, std::vector<Real> &d, std::vector<Real> &du,
213 const std::vector<Real> &r,
const bool transpose =
false) {
214 u.assign(r.begin(),r.end());
216 Teuchos::LAPACK<int,Real> lp;
217 std::vector<Real> du2(this->nx_-2,0.0);
218 std::vector<int> ipiv(this->nx_,0);
222 lp.GTTRF(this->nx_,&dl[0],&d[0],&du[0],&du2[0],&ipiv[0],&info);
227 lp.GTTRS(trans,this->nx_,nhrs,&dl[0],&d[0],&du[0],&du2[0],&ipiv[0],&u[0],ldb,&info);
230 void run_newton(std::vector<Real> &u,
const std::vector<Real> &z,
const std::vector<Real> ¶m) {
232 std::vector<Real> r(u.size(),0.0);
239 std::vector<Real> d(this->nx_,0.0);
240 std::vector<Real> dl(this->nx_-1,0.0);
241 std::vector<Real> du(this->nx_-1,0.0);
243 Real alpha = 1.0, tmp = 0.0;
244 std::vector<Real> s(this->nx_,0.0);
245 std::vector<Real> utmp(this->nx_,0.0);
246 for (
int i=0; i<maxit; i++) {
255 utmp.assign(u.begin(),u.end());
256 this->
update(utmp,s,-alpha);
259 while ( rnorm > (1.0-1.e-4*alpha)*tmp && alpha > std::sqrt(
ROL::ROL_EPSILON) ) {
261 utmp.assign(u.begin(),u.end());
262 this->
update(utmp,s,-alpha);
267 u.assign(utmp.begin(),utmp.end());
280 dx_ = 1.0/((Real)nx+1.0);
283 void solve_state(std::vector<Real> &u,
const std::vector<Real> &z,
const std::vector<Real> ¶m) {
285 u.resize(this->nx_,1.0);
289 void solve_adjoint(std::vector<Real> &p,
const std::vector<Real> &u,
const std::vector<Real> ¶m) {
294 std::vector<Real> d(this->nx_,0.0);
295 std::vector<Real> du(this->nx_-1,0.0);
296 std::vector<Real> dl(this->nx_-1,0.0);
299 std::vector<Real> r(this->nx_,0.0);
300 std::vector<Real> diff(this->nx_,0.0);
301 for (
int i=0; i<this->
nx_; i++) {
310 const std::vector<Real> &z,
const std::vector<Real> ¶m) {
315 std::vector<Real> d(this->nx_,0.0);
316 std::vector<Real> dl(this->nx_-1,0.0);
317 std::vector<Real> du(this->nx_-1,0.0);
320 std::vector<Real> r(this->nx_,0.0);
321 for (
int i=0; i<this->
nx_; i++) {
322 r[i] = this->dx_/6.0*(z[i]+4.0*z[i+1]+z[i+2]);
329 const std::vector<Real> &p,
const std::vector<Real> &v,
330 const std::vector<Real> &z,
const std::vector<Real> ¶m) {
335 std::vector<Real> d(this->nx_,0.0);
336 std::vector<Real> dl(this->nx_-1,0.0);
337 std::vector<Real> du(this->nx_-1,0.0);
340 std::vector<Real> r(this->nx_,0.0);
349 Teuchos::RCP<const std::vector<Real> > zp =
352 std::vector<Real> param(4,0.0);
356 Real val = this->alpha_*0.5*this->
dot(*zp,*zp);
357 Real res = 0.0, res1 = 0.0, res2 = 0.0, res3 = 0.0;
358 for (
int i=0; i<this->
nx_; i++) {
362 res = this->dx_/6.0*(4.0*res1 + res2)*res1;
364 else if ( i == this->nx_-1 ) {
367 res = this->dx_/6.0*(res1 + 4.0*res2)*res2;
373 res = this->dx_/6.0*(res1 + 4.0*res2 + res3)*res2;
381 Teuchos::RCP<const std::vector<Real> > zp =
383 Teuchos::RCP<std::vector<Real> > gp =
384 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<
ROL::StdVector<Real> >(g)).getVector());
386 std::vector<Real> param(4,0.0);
393 for (
int i=0; i<this->nx_+2; i++) {
395 (*gp)[i] = this->alpha_*this->dx_/6.0*(2.0*(*zp)[i]+(*zp)[i+1]);
396 (*gp)[i] -= this->dx_/6.0*(p[i]);
398 else if (i==this->nx_+1) {
399 (*gp)[i] = this->alpha_*this->dx_/6.0*(2.0*(*zp)[i]+(*zp)[i-1]);
400 (*gp)[i] -= this->dx_/6.0*(p[i-2]);
403 (*gp)[i] = this->alpha_*this->dx_/6.0*((*zp)[i-1]+4.0*(*zp)[i]+(*zp)[i+1]);
405 (*gp)[i] -= this->dx_/6.0*(4.0*p[i-1]+p[i]);
407 else if (i==this->nx_) {
408 (*gp)[i] -= this->dx_/6.0*(4.0*p[i-1]+p[i-2]);
411 (*gp)[i] -= this->dx_/6.0*(p[i-2]+4.0*p[i-1]+p[i]);
418 Teuchos::RCP<const std::vector<Real> > zp =
420 Teuchos::RCP<const std::vector<Real> > vp =
422 Teuchos::RCP<std::vector<Real> > hvp =
423 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<
ROL::StdVector<Real> >(hv)).getVector());
425 std::vector<Real> param(4,0.0);
438 for (
int i=0; i<this->nx_+2; i++) {
440 (*hvp)[i] = this->alpha_*this->dx_/6.0*(2.0*(*vp)[i]+(*vp)[i+1]);
441 (*hvp)[i] -= this->dx_/6.0*(q[i]);
443 else if (i==this->nx_+1) {
444 (*hvp)[i] = this->alpha_*this->dx_/6.0*(2.0*(*vp)[i]+(*vp)[i-1]);
445 (*hvp)[i] -= this->dx_/6.0*(q[i-2]);
448 (*hvp)[i] = this->alpha_*this->dx_/6.0*((*vp)[i-1]+4.0*(*vp)[i]+(*vp)[i+1]);
450 (*hvp)[i] -= this->dx_/6.0*(4.0*q[i-1]+q[i]);
452 else if (i==this->nx_) {
453 (*hvp)[i] -= this->dx_/6.0*(4.0*q[i-1]+q[i-2]);
456 (*hvp)[i] -= this->dx_/6.0*(q[i-2]+4.0*q[i-1]+q[i]);
472 x_lo_.resize(dim_,0.0);
473 x_up_.resize(dim_,1.0);
476 Teuchos::RCP<const std::vector<Real> > ex =
480 for (
int i = 0; i < this->dim_; i++ ) {
481 if ( (*ex)[i] >= this->x_lo_[i] && (*ex)[i] <= this->x_up_[i] ) { cnt *= 1; }
484 if ( cnt == 0 ) { val =
false; }
488 Teuchos::RCP<std::vector<Real> > ex =
489 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<
ROL::StdVector<Real> >(x)).getVector());
490 for (
int i = 0; i < this->dim_; i++ ) {
491 (*ex)[i] = std::max(this->x_lo_[i],std::min(this->x_up_[i],(*ex)[i]));
495 Teuchos::RCP<const std::vector<Real> > ex =
497 Teuchos::RCP<std::vector<Real> > ev =
498 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<
ROL::StdVector<Real> >(v)).getVector());
499 Real epsn = std::min(eps,this->min_diff_);
500 for (
int i = 0; i < this->dim_; i++ ) {
501 if ( ((*ex)[i] <= this->x_lo_[i]+epsn) ) {
507 Teuchos::RCP<const std::vector<Real> > ex =
509 Teuchos::RCP<std::vector<Real> > ev =
510 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<
ROL::StdVector<Real> >(v)).getVector());
511 Real epsn = std::min(eps,this->min_diff_);
512 for (
int i = 0; i < this->dim_; i++ ) {
513 if ( ((*ex)[i] >= this->x_up_[i]-epsn) ) {
519 Teuchos::RCP<const std::vector<Real> > ex =
521 Teuchos::RCP<const std::vector<Real> > eg =
523 Teuchos::RCP<std::vector<Real> > ev =
524 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<
ROL::StdVector<Real> >(v)).getVector());
525 Real epsn = std::min(eps,this->min_diff_);
526 for (
int i = 0; i < this->dim_; i++ ) {
527 if ( ((*ex)[i] <= this->x_lo_[i]+epsn && (*eg)[i] > 0.0) ){
533 Teuchos::RCP<const std::vector<Real> > ex =
535 Teuchos::RCP<const std::vector<Real> > eg =
537 Teuchos::RCP<std::vector<Real> > ev =
538 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<
ROL::StdVector<Real> >(v)).getVector());
539 Real epsn = std::min(eps,this->min_diff_);
540 for (
int i = 0; i < this->dim_; i++ ) {
541 if ( ((*ex)[i] >= this->x_up_[i]-epsn && (*eg)[i] < 0.0) ) {
547 Teuchos::RCP<std::vector<Real> > us = Teuchos::rcp(
new std::vector<Real>(this->dim_,0.0) );
548 us->assign(this->x_up_.begin(),this->x_up_.end());
553 Teuchos::RCP<std::vector<Real> > ls = Teuchos::rcp(
new std::vector<Real>(this->dim_,0.0) );
554 ls->assign(this->x_lo_.begin(),this->x_lo_.end());
Provides the interface to evaluate objective functions.
void pruneUpperActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &x, Real eps)
Set variables to zero if they correspond to the upper -active set.
Real evaluate_target(Real x)
void solve_adjoint_sensitivity(std::vector< Real > &q, const std::vector< Real > &u, const std::vector< Real > &p, const std::vector< Real > &v, const std::vector< Real > &z, const std::vector< Real > ¶m)
void solve_state(std::vector< Real > &u, const std::vector< Real > &z, const std::vector< Real > ¶m)
void run_newton(std::vector< Real > &u, const std::vector< Real > &z, const std::vector< Real > ¶m)
void pruneUpperActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &g, const ROL::Vector< Real > &x, Real eps)
Set variables to zero if they correspond to the upper -binding set.
Contains definitions of custom data types in ROL.
void solve_state_sensitivity(std::vector< Real > &v, const std::vector< Real > &u, const std::vector< Real > &z, const std::vector< Real > ¶m)
void compute_residual(std::vector< Real > &r, const std::vector< Real > &u, const std::vector< Real > &z, const std::vector< Real > ¶m)
Defines the linear algebra or vector space interface.
void pruneLowerActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &g, const ROL::Vector< Real > &x, Real eps)
Set variables to zero if they correspond to the lower -binding set.
void hessVec(ROL::Vector< Real > &hv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &z, Real &tol)
Apply Hessian approximation to vector.
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)
std::vector< Real > x_lo_
Real dot(const std::vector< Real > &x, const std::vector< Real > &y)
Objective_BurgersControl(Real alpha=1.e-4, int nx=128)
void pruneLowerActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &x, Real eps)
Set variables to zero if they correspond to the lower -active set.
void compute_pde_jacobian(std::vector< Real > &dl, std::vector< Real > &d, std::vector< Real > &du, const std::vector< Real > &u, const std::vector< Real > ¶m)
void update(std::vector< Real > &u, const std::vector< Real > &s, const Real alpha=1.0)
std::vector< Real > x_up_
void setVectorToLowerBound(ROL::Vector< Real > &l)
Set the input vector to the lower bound.
Provides the interface to apply upper and lower bound constraints.
void add_pde_hessian(std::vector< Real > &r, const std::vector< Real > &u, const std::vector< Real > &p, const std::vector< Real > &s, Real alpha=1.0)
void solve_adjoint(std::vector< Real > &p, const std::vector< Real > &u, const std::vector< Real > ¶m)
void project(ROL::Vector< Real > &x)
Project optimization variables onto the bounds.
Real value(const ROL::Vector< Real > &z, Real &tol)
Compute value.
void apply_mass(std::vector< Real > &Mu, const std::vector< Real > &u)
virtual void set(const Vector &x)
Set where .
void gradient(ROL::Vector< Real > &g, const ROL::Vector< Real > &z, Real &tol)
Compute gradient.
bool isFeasible(const ROL::Vector< Real > &x)
Check if the vector, v, is feasible.
BoundConstraint_BurgersControl(int dim)
void setVectorToUpperBound(ROL::Vector< Real > &u)
Set the input vector to the upper bound.
Real compute_norm(const std::vector< Real > &r)
void scale(std::vector< Real > &u, const Real alpha=0.0)
static const double ROL_EPSILON
Platform-dependent machine epsilon.