ROL
burgers-control/example_01.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ************************************************************************
3 //
4 // Rapid Optimization Library (ROL) Package
5 // Copyright (2014) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact lead developers:
38 // Drew Kouri (dpkouri@sandia.gov) and
39 // Denis Ridzal (dridzal@sandia.gov)
40 //
41 // ************************************************************************
42 // @HEADER
43 
49 #include "ROL_Algorithm.hpp"
50 #include "ROL_Types.hpp"
51 #include "Teuchos_oblackholestream.hpp"
52 #include "Teuchos_GlobalMPISession.hpp"
53 #include "Teuchos_XMLParameterListHelpers.hpp"
54 #include "Teuchos_LAPACK.hpp"
55 
56 #include <iostream>
57 #include <algorithm>
58 
59 #include "ROL_StdVector.hpp"
60 #include "ROL_Objective.hpp"
61 #include "ROL_BoundConstraint.hpp"
62 
63 template<class Real>
65 private:
66  Real alpha_; // Penalty Parameter
67 
68  int nx_; // Number of interior nodes
69  Real dx_; // Mesh spacing (i.e. 1/(nx+1))
70 
71 /***************************************************************/
72 /********** BEGIN PRIVATE MEMBER FUNCTION DECLARATION **********/
73 /***************************************************************/
74  Real evaluate_target(Real x) {
75  Real val = 0.0;
76  int example = 2;
77  switch (example) {
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;
82  }
83  return val;
84  }
85 
86  Real compute_norm(const std::vector<Real> &r) {
87  return std::sqrt(this->dot(r,r));
88  }
89 
90  Real dot(const std::vector<Real> &x, const std::vector<Real> &y) {
91  Real ip = 0.0;
92  Real c = (((int)x.size()==this->nx_) ? 4.0 : 2.0);
93  for (unsigned i=0; i<x.size(); i++) {
94  if ( i == 0 ) {
95  ip += this->dx_/6.0*(c*x[i] + x[i+1])*y[i];
96  }
97  else if ( i == x.size()-1 ) {
98  ip += this->dx_/6.0*(x[i-1] + c*x[i])*y[i];
99  }
100  else {
101  ip += this->dx_/6.0*(x[i-1] + 4.0*x[i] + x[i+1])*y[i];
102  }
103  }
104  return ip;
105  }
106 
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++) {
109  u[i] += alpha*s[i];
110  }
111  }
112 
113  void scale(std::vector<Real> &u, const Real alpha=0.0) {
114  for (unsigned i=0; i<u.size(); i++) {
115  u[i] *= alpha;
116  }
117  }
118 
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++) {
123  if ( i == 0 ) {
124  Mu[i] = this->dx_/6.0*(c*u[i] + u[i+1]);
125  }
126  else if ( i == u.size()-1 ) {
127  Mu[i] = this->dx_/6.0*(u[i-1] + c*u[i]);
128  }
129  else {
130  Mu[i] = this->dx_/6.0*(u[i-1] + 4.0*u[i] + u[i+1]);
131  }
132  }
133  }
134 
135  void compute_residual(std::vector<Real> &r, const std::vector<Real> &u,
136  const std::vector<Real> &z, const std::vector<Real> &param) {
137  r.clear();
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++) {
144  // Contribution from stiffness term
145  if ( i==0 ) {
146  r[i] = nu/this->dx_*(2.0*u[i]-u[i+1]);
147  }
148  else if (i==this->nx_-1) {
149  r[i] = nu/this->dx_*(2.0*u[i]-u[i-1]);
150  }
151  else {
152  r[i] = nu/this->dx_*(2.0*u[i]-u[i-1]-u[i+1]);
153  }
154  // Contribution from nonlinear term
155  if (i<this->nx_-1){
156  r[i] += u[i+1]*(u[i]+u[i+1])/6.0;
157  }
158  if (i>0) {
159  r[i] -= u[i-1]*(u[i-1]+u[i])/6.0;
160  }
161  // Contribution from control
162  r[i] -= this->dx_/6.0*(z[i]+4.0*z[i+1]+z[i+2]);
163  // Contribution from right-hand side
164  r[i] -= this->dx_*f;
165  }
166  // Contribution from Dirichlet boundary terms
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_;
169  }
170 
171  void compute_pde_jacobian(std::vector<Real> &dl, std::vector<Real> &d, std::vector<Real> &du,
172  const std::vector<Real> &u, const std::vector<Real> &param) {
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;
176  // Get Diagonal and Off-Diagonal Entries of linear PDE Jacobian
177  d.clear();
178  d.resize(this->nx_,nu*2.0/this->dx_);
179  dl.clear();
180  dl.resize(this->nx_-1,-nu/this->dx_);
181  du.clear();
182  du.resize(this->nx_-1,-nu/this->dx_);
183  // Contribution from nonlinearity
184  for (int i=0; i<this->nx_; i++) {
185  if (i<this->nx_-1) {
186  dl[i] += (-2.0*u[i]-u[i+1])/6.0;
187  d[i] += u[i+1]/6.0;
188  }
189  if (i>0) {
190  d[i] += -u[i-1]/6.0;
191  du[i-1] += (u[i-1]+2.0*u[i])/6.0;
192  }
193  }
194  // Contribution from Dirichlet boundary conditions
195  d[0] -= u0/6.0;
196  d[this->nx_-1] += u1/6.0;
197  }
198 
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++) {
202  // Contribution from nonlinear term
203  if (i<this->nx_-1){
204  r[i] += alpha*(p[i]*s[i+1] - p[i+1]*(2.0*s[i]+s[i+1]))/6.0;
205  }
206  if (i>0) {
207  r[i] += alpha*(p[i-1]*(s[i-1]+2.0*s[i]) - p[i]*s[i-1])/6.0;
208  }
209  }
210  }
211 
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());
215  // Perform LDL factorization
216  Teuchos::LAPACK<int,Real> lp;
217  std::vector<Real> du2(this->nx_-2,0.0);
218  std::vector<int> ipiv(this->nx_,0);
219  int info;
220  int ldb = this->nx_;
221  int nhrs = 1;
222  lp.GTTRF(this->nx_,&dl[0],&d[0],&du[0],&du2[0],&ipiv[0],&info);
223  char trans = 'N';
224  if ( transpose ) {
225  trans = 'T';
226  }
227  lp.GTTRS(trans,this->nx_,nhrs,&dl[0],&d[0],&du[0],&du2[0],&ipiv[0],&u[0],ldb,&info);
228  }
229 
230  void run_newton(std::vector<Real> &u, const std::vector<Real> &z, const std::vector<Real> &param) {
231  // Compute residual and residual norm
232  std::vector<Real> r(u.size(),0.0);
233  this->compute_residual(r,u,z,param);
234  Real rnorm = this->compute_norm(r);
235  // Define tolerances
236  Real tol = 1.e2*ROL::ROL_EPSILON;
237  Real maxit = 500;
238  // Initialize Jacobian storage
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);
242  // Iterate Newton's method
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++) {
247  //std::cout << i << " " << rnorm << "\n";
248  // Get Jacobian
249  this->compute_pde_jacobian(dl,d,du,u,param);
250  // Solve Newton system
251  this->linear_solve(s,dl,d,du,r);
252  // Perform line search
253  tmp = rnorm;
254  alpha = 1.0;
255  utmp.assign(u.begin(),u.end());
256  this->update(utmp,s,-alpha);
257  this->compute_residual(r,utmp,z,param);
258  rnorm = this->compute_norm(r);
259  while ( rnorm > (1.0-1.e-4*alpha)*tmp && alpha > std::sqrt(ROL::ROL_EPSILON) ) {
260  alpha /= 2.0;
261  utmp.assign(u.begin(),u.end());
262  this->update(utmp,s,-alpha);
263  this->compute_residual(r,utmp,z,param);
264  rnorm = this->compute_norm(r);
265  }
266  // Update iterate
267  u.assign(utmp.begin(),utmp.end());
268  if ( rnorm < tol ) {
269  break;
270  }
271  }
272  }
273 /*************************************************************/
274 /********** END PRIVATE MEMBER FUNCTION DECLARATION **********/
275 /*************************************************************/
276 
277 public:
278 
279  Objective_BurgersControl(Real alpha = 1.e-4, int nx = 128) : alpha_(alpha), nx_(nx) {
280  dx_ = 1.0/((Real)nx+1.0);
281  }
282 
283  void solve_state(std::vector<Real> &u, const std::vector<Real> &z, const std::vector<Real> &param) {
284  u.clear();
285  u.resize(this->nx_,1.0);
286  this->run_newton(u,z,param);
287  }
288 
289  void solve_adjoint(std::vector<Real> &p, const std::vector<Real> &u, const std::vector<Real> &param) {
290  // Initialize State Storage
291  p.clear();
292  p.resize(this->nx_);
293  // Get PDE Jacobian
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);
297  this->compute_pde_jacobian(dl,d,du,u,param);
298  // Get Right Hand Side
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++) {
302  diff[i] = -(u[i]-this->evaluate_target((Real)(i+1)*this->dx_));
303  }
304  this->apply_mass(r,diff);
305  // Solve solve adjoint system at current time step
306  this->linear_solve(p,dl,d,du,r,true);
307  }
308 
309  void solve_state_sensitivity(std::vector<Real> &v, const std::vector<Real> &u,
310  const std::vector<Real> &z, const std::vector<Real> &param) {
311  // Initialize State Storage
312  v.clear();
313  v.resize(this->nx_);
314  // Get PDE Jacobian
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);
318  this->compute_pde_jacobian(dl,d,du,u,param);
319  // Get Right Hand Side
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]);
323  }
324  // Solve solve state sensitivity system at current time step
325  this->linear_solve(v,dl,d,du,r);
326  }
327 
328  void solve_adjoint_sensitivity(std::vector<Real> &q, const std::vector<Real> &u,
329  const std::vector<Real> &p, const std::vector<Real> &v,
330  const std::vector<Real> &z, const std::vector<Real> &param) {
331  // Initialize State Storage
332  q.clear();
333  q.resize(this->nx_);
334  // Get PDE Jacobian
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);
338  this->compute_pde_jacobian(dl,d,du,u,param);
339  // Get Right Hand Side
340  std::vector<Real> r(this->nx_,0.0);
341  this->apply_mass(r,v);
342  this->scale(r,-1.0);
343  this->add_pde_hessian(r,u,p,v,-1.0);
344  // Solve Tridiagonal System Using LAPACK's SPD Tridiagonal Solver
345  this->linear_solve(q,dl,d,du,r,true);
346  }
347 
348  Real value( const ROL::Vector<Real> &z, Real &tol ) {
349  Teuchos::RCP<const std::vector<Real> > zp =
350  (Teuchos::dyn_cast<ROL::StdVector<Real> >(const_cast<ROL::Vector<Real> &>(z))).getVector();
351  // SOLVE STATE EQUATION
352  std::vector<Real> param(4,0.0);
353  std::vector<Real> u;
354  this->solve_state(u,*zp,param);
355  // COMPUTE RESIDUAL
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++) {
359  if ( i == 0 ) {
360  res1 = u[i]-evaluate_target((Real)(i+1)*this->dx_);
361  res2 = u[i+1]-evaluate_target((Real)(i+2)*this->dx_);
362  res = this->dx_/6.0*(4.0*res1 + res2)*res1;
363  }
364  else if ( i == this->nx_-1 ) {
365  res1 = u[i-1]-evaluate_target((Real)i*this->dx_);
366  res2 = u[i]-evaluate_target((Real)(i+1)*this->dx_);
367  res = this->dx_/6.0*(res1 + 4.0*res2)*res2;
368  }
369  else {
370  res1 = u[i-1]-evaluate_target((Real)i*this->dx_);
371  res2 = u[i]-evaluate_target((Real)(i+1)*this->dx_);
372  res3 = u[i+1]-evaluate_target((Real)(i+2)*this->dx_);
373  res = this->dx_/6.0*(res1 + 4.0*res2 + res3)*res2;
374  }
375  val += 0.5*res;
376  }
377  return val;
378  }
379 
380  void gradient( ROL::Vector<Real> &g, const ROL::Vector<Real> &z, Real &tol ) {
381  Teuchos::RCP<const std::vector<Real> > zp =
382  (Teuchos::dyn_cast<ROL::StdVector<Real> >(const_cast<ROL::Vector<Real> &>(z))).getVector();
383  Teuchos::RCP<std::vector<Real> > gp =
384  Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<ROL::StdVector<Real> >(g)).getVector());
385  // SOLVE STATE EQUATION
386  std::vector<Real> param(4,0.0);
387  std::vector<Real> u;
388  this->solve_state(u,*zp,param);
389  // SOLVE ADJOINT EQUATION
390  std::vector<Real> p;
391  this->solve_adjoint(p,u,param);
392  // COMPUTE GRADIENT
393  for (int i=0; i<this->nx_+2; i++) {
394  if (i==0) {
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]);
397  }
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]);
401  }
402  else {
403  (*gp)[i] = this->alpha_*this->dx_/6.0*((*zp)[i-1]+4.0*(*zp)[i]+(*zp)[i+1]);
404  if (i==1) {
405  (*gp)[i] -= this->dx_/6.0*(4.0*p[i-1]+p[i]);
406  }
407  else if (i==this->nx_) {
408  (*gp)[i] -= this->dx_/6.0*(4.0*p[i-1]+p[i-2]);
409  }
410  else {
411  (*gp)[i] -= this->dx_/6.0*(p[i-2]+4.0*p[i-1]+p[i]);
412  }
413  }
414  }
415  }
416 
417  void hessVec( ROL::Vector<Real> &hv, const ROL::Vector<Real> &v, const ROL::Vector<Real> &z, Real &tol ) {
418  Teuchos::RCP<const std::vector<Real> > zp =
419  (Teuchos::dyn_cast<ROL::StdVector<Real> >(const_cast<ROL::Vector<Real> &>(z))).getVector();
420  Teuchos::RCP<const std::vector<Real> > vp =
421  (Teuchos::dyn_cast<ROL::StdVector<Real> >(const_cast<ROL::Vector<Real> &>(v))).getVector();
422  Teuchos::RCP<std::vector<Real> > hvp =
423  Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<ROL::StdVector<Real> >(hv)).getVector());
424  // SOLVE STATE EQUATION
425  std::vector<Real> param(4,0.0);
426  std::vector<Real> u;
427  this->solve_state(u,*zp,param);
428  // SOLVE ADJOINT EQUATION
429  std::vector<Real> p;
430  this->solve_adjoint(p,u,param);
431  // SOLVE STATE SENSITIVITY EQUATION
432  std::vector<Real> s;
433  this->solve_state_sensitivity(s,u,*vp,param);
434  // SOLVE ADJOINT SENSITIVITY EQUATION
435  std::vector<Real> q;
436  this->solve_adjoint_sensitivity(q,u,p,s,*vp,param);
437  // COMPUTE HESSVEC
438  for (int i=0; i<this->nx_+2; i++) {
439  if (i==0) {
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]);
442  }
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]);
446  }
447  else {
448  (*hvp)[i] = this->alpha_*this->dx_/6.0*((*vp)[i-1]+4.0*(*vp)[i]+(*vp)[i+1]);
449  if (i==1) {
450  (*hvp)[i] -= this->dx_/6.0*(4.0*q[i-1]+q[i]);
451  }
452  else if (i==this->nx_) {
453  (*hvp)[i] -= this->dx_/6.0*(4.0*q[i-1]+q[i-2]);
454  }
455  else {
456  (*hvp)[i] -= this->dx_/6.0*(q[i-2]+4.0*q[i-1]+q[i]);
457  }
458  }
459  }
460  }
461 };
462 
463 template<class Real>
465 private:
466  int dim_;
467  std::vector<Real> x_lo_;
468  std::vector<Real> x_up_;
469  Real min_diff_;
470 public:
471  BoundConstraint_BurgersControl(int dim) : dim_(dim), min_diff_(0.5) {
472  x_lo_.resize(dim_,0.0);
473  x_up_.resize(dim_,1.0);
474  }
475  bool isFeasible( const ROL::Vector<Real> &x ) {
476  Teuchos::RCP<const std::vector<Real> > ex =
477  (Teuchos::dyn_cast<ROL::StdVector<Real> >(const_cast<ROL::Vector<Real> &>(x))).getVector();
478  bool val = true;
479  int cnt = 1;
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; }
482  else { cnt *= 0; }
483  }
484  if ( cnt == 0 ) { val = false; }
485  return val;
486  }
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]));
492  }
493  }
494  void pruneLowerActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &x, Real eps) {
495  Teuchos::RCP<const std::vector<Real> > ex =
496  (Teuchos::dyn_cast<ROL::StdVector<Real> >(const_cast<ROL::Vector<Real> &>(x))).getVector();
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) ) {
502  (*ev)[i] = 0.0;
503  }
504  }
505  }
506  void pruneUpperActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &x, Real eps) {
507  Teuchos::RCP<const std::vector<Real> > ex =
508  (Teuchos::dyn_cast<ROL::StdVector<Real> >(const_cast<ROL::Vector<Real> &>(x))).getVector();
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) ) {
514  (*ev)[i] = 0.0;
515  }
516  }
517  }
519  Teuchos::RCP<const std::vector<Real> > ex =
520  (Teuchos::dyn_cast<ROL::StdVector<Real> >(const_cast<ROL::Vector<Real> &>(x))).getVector();
521  Teuchos::RCP<const std::vector<Real> > eg =
522  (Teuchos::dyn_cast<ROL::StdVector<Real> >(const_cast<ROL::Vector<Real> &>(g))).getVector();
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) ){
528  (*ev)[i] = 0.0;
529  }
530  }
531  }
533  Teuchos::RCP<const std::vector<Real> > ex =
534  (Teuchos::dyn_cast<ROL::StdVector<Real> >(const_cast<ROL::Vector<Real> &>(x))).getVector();
535  Teuchos::RCP<const std::vector<Real> > eg =
536  (Teuchos::dyn_cast<ROL::StdVector<Real> >(const_cast<ROL::Vector<Real> &>(g))).getVector();
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) ) {
542  (*ev)[i] = 0.0;
543  }
544  }
545  }
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());
549  Teuchos::RCP<ROL::Vector<Real> > up = Teuchos::rcp( new ROL::StdVector<Real>(us) );
550  u.set(*up);
551  }
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());
555  Teuchos::RCP<ROL::Vector<Real> > lp = Teuchos::rcp( new ROL::StdVector<Real>(ls) );
556  l.set(*lp);
557  }
558 };
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.
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 > &param)
void solve_state(std::vector< Real > &u, const std::vector< Real > &z, const std::vector< Real > &param)
void run_newton(std::vector< Real > &u, const std::vector< Real > &z, const std::vector< Real > &param)
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 > &param)
void compute_residual(std::vector< Real > &r, const std::vector< Real > &u, const std::vector< Real > &z, const std::vector< Real > &param)
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:74
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)
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 > &param)
void update(std::vector< Real > &u, const std::vector< Real > &s, const Real alpha=1.0)
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 > &param)
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 .
Definition: ROL_Vector.hpp:196
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.
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.
Definition: ROL_Types.hpp:118