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 "ROL_Stream.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_Bounds.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 
108 
109  void update(std::vector<Real> &u, const std::vector<Real> &s, const Real alpha=1.0) {
110  for (unsigned i=0; i<u.size(); i++) {
111  u[i] += alpha*s[i];
112  }
113  }
114 
115  void scale(std::vector<Real> &u, const Real alpha=0.0) {
116  for (unsigned i=0; i<u.size(); i++) {
117  u[i] *= alpha;
118  }
119  }
120 
121  void apply_mass(std::vector<Real> &Mu, const std::vector<Real> &u ) {
122  Mu.resize(u.size(),0.0);
123  Real c = (((int)u.size()==this->nx_) ? 4.0 : 2.0);
124  for (unsigned i=0; i<u.size(); i++) {
125  if ( i == 0 ) {
126  Mu[i] = this->dx_/6.0*(c*u[i] + u[i+1]);
127  }
128  else if ( i == u.size()-1 ) {
129  Mu[i] = this->dx_/6.0*(u[i-1] + c*u[i]);
130  }
131  else {
132  Mu[i] = this->dx_/6.0*(u[i-1] + 4.0*u[i] + u[i+1]);
133  }
134  }
135  }
136 
137  void compute_residual(std::vector<Real> &r, const std::vector<Real> &u,
138  const std::vector<Real> &z, const std::vector<Real> &param) {
139  r.clear();
140  r.resize(this->nx_,0.0);
141  Real nu = std::pow(10.0,param[0]-2.0);
142  Real f = param[1]/100.0;
143  Real u0 = 1.0+param[2]/1000.0;
144  Real u1 = param[3]/1000.0;
145  for (int i=0; i<this->nx_; i++) {
146  // Contribution from stiffness term
147  if ( i==0 ) {
148  r[i] = nu/this->dx_*(2.0*u[i]-u[i+1]);
149  }
150  else if (i==this->nx_-1) {
151  r[i] = nu/this->dx_*(2.0*u[i]-u[i-1]);
152  }
153  else {
154  r[i] = nu/this->dx_*(2.0*u[i]-u[i-1]-u[i+1]);
155  }
156  // Contribution from nonlinear term
157  if (i<this->nx_-1){
158  r[i] += u[i+1]*(u[i]+u[i+1])/6.0;
159  }
160  if (i>0) {
161  r[i] -= u[i-1]*(u[i-1]+u[i])/6.0;
162  }
163  // Contribution from control
164  r[i] -= this->dx_/6.0*(z[i]+4.0*z[i+1]+z[i+2]);
165  // Contribution from right-hand side
166  r[i] -= this->dx_*f;
167  }
168  // Contribution from Dirichlet boundary terms
169  r[0] -= u0*u[ 0]/6.0 + u0*u0/6.0 + nu*u0/this->dx_;
170  r[this->nx_-1] += u1*u[this->nx_-1]/6.0 + u1*u1/6.0 - nu*u1/this->dx_;
171  }
172 
173  void compute_pde_jacobian(std::vector<Real> &dl, std::vector<Real> &d, std::vector<Real> &du,
174  const std::vector<Real> &u, const std::vector<Real> &param) {
175  Real nu = std::pow(10.0,param[0]-2.0);
176  Real u0 = 1.0+param[2]/1000.0;
177  Real u1 = param[3]/1000.0;
178  // Get Diagonal and Off-Diagonal Entries of linear PDE Jacobian
179  d.clear();
180  d.resize(this->nx_,nu*2.0/this->dx_);
181  dl.clear();
182  dl.resize(this->nx_-1,-nu/this->dx_);
183  du.clear();
184  du.resize(this->nx_-1,-nu/this->dx_);
185  // Contribution from nonlinearity
186  for (int i=0; i<this->nx_; i++) {
187  if (i<this->nx_-1) {
188  dl[i] += (-2.0*u[i]-u[i+1])/6.0;
189  d[i] += u[i+1]/6.0;
190  }
191  if (i>0) {
192  d[i] += -u[i-1]/6.0;
193  du[i-1] += (u[i-1]+2.0*u[i])/6.0;
194  }
195  }
196  // Contribution from Dirichlet boundary conditions
197  d[0] -= u0/6.0;
198  d[this->nx_-1] += u1/6.0;
199  }
200 
201  void add_pde_hessian(std::vector<Real> &r, const std::vector<Real> &u, const std::vector<Real> &p,
202  const std::vector<Real> &s, Real alpha = 1.0) {
203  for (int i=0; i<this->nx_; i++) {
204  // Contribution from nonlinear term
205  if (i<this->nx_-1){
206  r[i] += alpha*(p[i]*s[i+1] - p[i+1]*(2.0*s[i]+s[i+1]))/6.0;
207  }
208  if (i>0) {
209  r[i] += alpha*(p[i-1]*(s[i-1]+2.0*s[i]) - p[i]*s[i-1])/6.0;
210  }
211  }
212  }
213 
214  void linear_solve(std::vector<Real> &u, std::vector<Real> &dl, std::vector<Real> &d, std::vector<Real> &du,
215  const std::vector<Real> &r, const bool transpose = false) {
216  u.assign(r.begin(),r.end());
217  // Perform LDL factorization
218  Teuchos::LAPACK<int,Real> lp;
219  std::vector<Real> du2(this->nx_-2,0.0);
220  std::vector<int> ipiv(this->nx_,0);
221  int info;
222  int ldb = this->nx_;
223  int nhrs = 1;
224  lp.GTTRF(this->nx_,&dl[0],&d[0],&du[0],&du2[0],&ipiv[0],&info);
225  char trans = 'N';
226  if ( transpose ) {
227  trans = 'T';
228  }
229  lp.GTTRS(trans,this->nx_,nhrs,&dl[0],&d[0],&du[0],&du2[0],&ipiv[0],&u[0],ldb,&info);
230  }
231 
232  void run_newton(std::vector<Real> &u, const std::vector<Real> &z, const std::vector<Real> &param) {
233  // Compute residual and residual norm
234  std::vector<Real> r(u.size(),0.0);
235  this->compute_residual(r,u,z,param);
236  Real rnorm = this->compute_norm(r);
237  // Define tolerances
238  Real tol = 1.e2*ROL::ROL_EPSILON<Real>();
239  Real maxit = 500;
240  // Initialize Jacobian storage
241  std::vector<Real> d(this->nx_,0.0);
242  std::vector<Real> dl(this->nx_-1,0.0);
243  std::vector<Real> du(this->nx_-1,0.0);
244  // Iterate Newton's method
245  Real alpha = 1.0, tmp = 0.0;
246  std::vector<Real> s(this->nx_,0.0);
247  std::vector<Real> utmp(this->nx_,0.0);
248  for (int i=0; i<maxit; i++) {
249  //std::cout << i << " " << rnorm << "\n";
250  // Get Jacobian
251  this->compute_pde_jacobian(dl,d,du,u,param);
252  // Solve Newton system
253  this->linear_solve(s,dl,d,du,r);
254  // Perform line search
255  tmp = rnorm;
256  alpha = 1.0;
257  utmp.assign(u.begin(),u.end());
258  this->update(utmp,s,-alpha);
259  this->compute_residual(r,utmp,z,param);
260  rnorm = this->compute_norm(r);
261  while ( rnorm > (1.0-1.e-4*alpha)*tmp && alpha > std::sqrt(ROL::ROL_EPSILON<Real>()) ) {
262  alpha /= 2.0;
263  utmp.assign(u.begin(),u.end());
264  this->update(utmp,s,-alpha);
265  this->compute_residual(r,utmp,z,param);
266  rnorm = this->compute_norm(r);
267  }
268  // Update iterate
269  u.assign(utmp.begin(),utmp.end());
270  if ( rnorm < tol ) {
271  break;
272  }
273  }
274  }
275 /*************************************************************/
276 /********** END PRIVATE MEMBER FUNCTION DECLARATION **********/
277 /*************************************************************/
278 
279 public:
280 
281  Objective_BurgersControl(Real alpha = 1.e-4, int nx = 128) : alpha_(alpha), nx_(nx) {
282  dx_ = 1.0/((Real)nx+1.0);
283  }
284 
285  void solve_state(std::vector<Real> &u, const std::vector<Real> &z, const std::vector<Real> &param) {
286  u.clear();
287  u.resize(this->nx_,1.0);
288  this->run_newton(u,z,param);
289  }
290 
291  void solve_adjoint(std::vector<Real> &p, const std::vector<Real> &u, const std::vector<Real> &param) {
292  // Initialize State Storage
293  p.clear();
294  p.resize(this->nx_);
295  // Get PDE Jacobian
296  std::vector<Real> d(this->nx_,0.0);
297  std::vector<Real> du(this->nx_-1,0.0);
298  std::vector<Real> dl(this->nx_-1,0.0);
299  this->compute_pde_jacobian(dl,d,du,u,param);
300  // Get Right Hand Side
301  std::vector<Real> r(this->nx_,0.0);
302  std::vector<Real> diff(this->nx_,0.0);
303  for (int i=0; i<this->nx_; i++) {
304  diff[i] = -(u[i]-this->evaluate_target((Real)(i+1)*this->dx_));
305  }
306  this->apply_mass(r,diff);
307  // Solve solve adjoint system at current time step
308  this->linear_solve(p,dl,d,du,r,true);
309  }
310 
311  void solve_state_sensitivity(std::vector<Real> &v, const std::vector<Real> &u,
312  const std::vector<Real> &z, const std::vector<Real> &param) {
313  // Initialize State Storage
314  v.clear();
315  v.resize(this->nx_);
316  // Get PDE Jacobian
317  std::vector<Real> d(this->nx_,0.0);
318  std::vector<Real> dl(this->nx_-1,0.0);
319  std::vector<Real> du(this->nx_-1,0.0);
320  this->compute_pde_jacobian(dl,d,du,u,param);
321  // Get Right Hand Side
322  std::vector<Real> r(this->nx_,0.0);
323  for (int i=0; i<this->nx_; i++) {
324  r[i] = this->dx_/6.0*(z[i]+4.0*z[i+1]+z[i+2]);
325  }
326  // Solve solve state sensitivity system at current time step
327  this->linear_solve(v,dl,d,du,r);
328  }
329 
330  void solve_adjoint_sensitivity(std::vector<Real> &q, const std::vector<Real> &u,
331  const std::vector<Real> &p, const std::vector<Real> &v,
332  const std::vector<Real> &z, const std::vector<Real> &param) {
333  // Initialize State Storage
334  q.clear();
335  q.resize(this->nx_);
336  // Get PDE Jacobian
337  std::vector<Real> d(this->nx_,0.0);
338  std::vector<Real> dl(this->nx_-1,0.0);
339  std::vector<Real> du(this->nx_-1,0.0);
340  this->compute_pde_jacobian(dl,d,du,u,param);
341  // Get Right Hand Side
342  std::vector<Real> r(this->nx_,0.0);
343  this->apply_mass(r,v);
344  this->scale(r,-1.0);
345  this->add_pde_hessian(r,u,p,v,-1.0);
346  // Solve Tridiagonal System Using LAPACK's SPD Tridiagonal Solver
347  this->linear_solve(q,dl,d,du,r,true);
348  }
349 
350  Real value( const ROL::Vector<Real> &z, Real &tol ) {
351  ROL::Ptr<const std::vector<Real> > zp =
352  dynamic_cast<const ROL::StdVector<Real>&>(z).getVector();
353  // SOLVE STATE EQUATION
354  std::vector<Real> param(4,0.0);
355  std::vector<Real> u;
356  this->solve_state(u,*zp,param);
357  // COMPUTE RESIDUAL
358  Real val = this->alpha_*0.5*this->dot(*zp,*zp);
359  Real res = 0.0, res1 = 0.0, res2 = 0.0, res3 = 0.0;
360  for (int i=0; i<this->nx_; i++) {
361  if ( i == 0 ) {
362  res1 = u[i]-evaluate_target((Real)(i+1)*this->dx_);
363  res2 = u[i+1]-evaluate_target((Real)(i+2)*this->dx_);
364  res = this->dx_/6.0*(4.0*res1 + res2)*res1;
365  }
366  else if ( i == this->nx_-1 ) {
367  res1 = u[i-1]-evaluate_target((Real)i*this->dx_);
368  res2 = u[i]-evaluate_target((Real)(i+1)*this->dx_);
369  res = this->dx_/6.0*(res1 + 4.0*res2)*res2;
370  }
371  else {
372  res1 = u[i-1]-evaluate_target((Real)i*this->dx_);
373  res2 = u[i]-evaluate_target((Real)(i+1)*this->dx_);
374  res3 = u[i+1]-evaluate_target((Real)(i+2)*this->dx_);
375  res = this->dx_/6.0*(res1 + 4.0*res2 + res3)*res2;
376  }
377  val += 0.5*res;
378  }
379  return val;
380  }
381 
382  void gradient( ROL::Vector<Real> &g, const ROL::Vector<Real> &z, Real &tol ) {
383  ROL::Ptr<const std::vector<Real> > zp =
384  dynamic_cast<const ROL::StdVector<Real>&>(z).getVector();
385  ROL::Ptr<std::vector<Real> > gp =
386  dynamic_cast<ROL::StdVector<Real>&>(g).getVector();
387  // SOLVE STATE EQUATION
388  std::vector<Real> param(4,0.0);
389  std::vector<Real> u;
390  this->solve_state(u,*zp,param);
391  // SOLVE ADJOINT EQUATION
392  std::vector<Real> p;
393  this->solve_adjoint(p,u,param);
394  // COMPUTE GRADIENT
395  for (int i=0; i<this->nx_+2; i++) {
396  if (i==0) {
397  (*gp)[i] = this->alpha_*this->dx_/6.0*(2.0*(*zp)[i]+(*zp)[i+1]);
398  (*gp)[i] -= this->dx_/6.0*(p[i]);
399  }
400  else if (i==this->nx_+1) {
401  (*gp)[i] = this->alpha_*this->dx_/6.0*(2.0*(*zp)[i]+(*zp)[i-1]);
402  (*gp)[i] -= this->dx_/6.0*(p[i-2]);
403  }
404  else {
405  (*gp)[i] = this->alpha_*this->dx_/6.0*((*zp)[i-1]+4.0*(*zp)[i]+(*zp)[i+1]);
406  if (i==1) {
407  (*gp)[i] -= this->dx_/6.0*(4.0*p[i-1]+p[i]);
408  }
409  else if (i==this->nx_) {
410  (*gp)[i] -= this->dx_/6.0*(4.0*p[i-1]+p[i-2]);
411  }
412  else {
413  (*gp)[i] -= this->dx_/6.0*(p[i-2]+4.0*p[i-1]+p[i]);
414  }
415  }
416  }
417  }
418 
419  void hessVec( ROL::Vector<Real> &hv, const ROL::Vector<Real> &v, const ROL::Vector<Real> &z, Real &tol ) {
420  ROL::Ptr<const std::vector<Real> > zp =
421  dynamic_cast<const ROL::StdVector<Real>&>(z).getVector();
422  ROL::Ptr<const std::vector<Real> > vp =
423  dynamic_cast<const ROL::StdVector<Real>&>(v).getVector();
424  ROL::Ptr<std::vector<Real> > hvp =
425  dynamic_cast<ROL::StdVector<Real>&>(hv).getVector();
426  // SOLVE STATE EQUATION
427  std::vector<Real> param(4,0.0);
428  std::vector<Real> u;
429  this->solve_state(u,*zp,param);
430  // SOLVE ADJOINT EQUATION
431  std::vector<Real> p;
432  this->solve_adjoint(p,u,param);
433  // SOLVE STATE SENSITIVITY EQUATION
434  std::vector<Real> s;
435  this->solve_state_sensitivity(s,u,*vp,param);
436  // SOLVE ADJOINT SENSITIVITY EQUATION
437  std::vector<Real> q;
438  this->solve_adjoint_sensitivity(q,u,p,s,*vp,param);
439  // COMPUTE HESSVEC
440  for (int i=0; i<this->nx_+2; i++) {
441  if (i==0) {
442  (*hvp)[i] = this->alpha_*this->dx_/6.0*(2.0*(*vp)[i]+(*vp)[i+1]);
443  (*hvp)[i] -= this->dx_/6.0*(q[i]);
444  }
445  else if (i==this->nx_+1) {
446  (*hvp)[i] = this->alpha_*this->dx_/6.0*(2.0*(*vp)[i]+(*vp)[i-1]);
447  (*hvp)[i] -= this->dx_/6.0*(q[i-2]);
448  }
449  else {
450  (*hvp)[i] = this->alpha_*this->dx_/6.0*((*vp)[i-1]+4.0*(*vp)[i]+(*vp)[i+1]);
451  if (i==1) {
452  (*hvp)[i] -= this->dx_/6.0*(4.0*q[i-1]+q[i]);
453  }
454  else if (i==this->nx_) {
455  (*hvp)[i] -= this->dx_/6.0*(4.0*q[i-1]+q[i-2]);
456  }
457  else {
458  (*hvp)[i] -= this->dx_/6.0*(q[i-2]+4.0*q[i-1]+q[i]);
459  }
460  }
461  }
462  }
463 };
464 
465 template<class Real>
467 private:
468  int dim_;
469  std::vector<Real> x_lo_;
470  std::vector<Real> x_up_;
471  Real min_diff_;
472 public:
474  x_lo_.resize(dim_,0.0);
475  x_up_.resize(dim_,1.0);
476  }
477  bool isFeasible( const ROL::Vector<Real> &x ) {
478  ROL::Ptr<const std::vector<Real> > ex =
479  dynamic_cast<const ROL::StdVector<Real>&>(x).getVector();
480  bool val = true;
481  int cnt = 1;
482  for ( int i = 0; i < this->dim_; i++ ) {
483  if ( (*ex)[i] >= this->x_lo_[i] && (*ex)[i] <= this->x_up_[i] ) { cnt *= 1; }
484  else { cnt *= 0; }
485  }
486  if ( cnt == 0 ) { val = false; }
487  return val;
488  }
490  ROL::Ptr<std::vector<Real> > ex =
491  dynamic_cast<ROL::StdVector<Real>&>(x).getVector();
492  for ( int i = 0; i < this->dim_; i++ ) {
493  (*ex)[i] = std::max(this->x_lo_[i],std::min(this->x_up_[i],(*ex)[i]));
494  }
495  }
496  void pruneLowerActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &x, Real eps) {
497  ROL::Ptr<const std::vector<Real> > ex =
498  dynamic_cast<const ROL::StdVector<Real>&>(x).getVector();
499  ROL::Ptr<std::vector<Real> > ev =
500  dynamic_cast<ROL::StdVector<Real>&>(v).getVector();
501  Real epsn = std::min(eps,this->min_diff_);
502  for ( int i = 0; i < this->dim_; i++ ) {
503  if ( ((*ex)[i] <= this->x_lo_[i]+epsn) ) {
504  (*ev)[i] = 0.0;
505  }
506  }
507  }
508  void pruneUpperActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &x, Real eps) {
509  ROL::Ptr<const std::vector<Real> > ex =
510  dynamic_cast<const ROL::StdVector<Real>&>(x).getVector();
511  ROL::Ptr<std::vector<Real> > ev =
512  dynamic_cast<ROL::StdVector<Real>&>(v).getVector();
513  Real epsn = std::min(eps,this->min_diff_);
514  for ( int i = 0; i < this->dim_; i++ ) {
515  if ( ((*ex)[i] >= this->x_up_[i]-epsn) ) {
516  (*ev)[i] = 0.0;
517  }
518  }
519  }
521  ROL::Ptr<const std::vector<Real> > ex =
522  dynamic_cast<ROL::StdVector<Real>&>(x).getVector();
523  ROL::Ptr<const std::vector<Real> > eg =
524  dynamic_cast<const ROL::StdVector<Real>&>(g).getVector();
525  ROL::Ptr<std::vector<Real> > ev =
526  dynamic_cast<ROL::StdVector<Real>&>(v).getVector();
527  Real epsn = std::min(eps,this->min_diff_);
528  for ( int i = 0; i < this->dim_; i++ ) {
529  if ( ((*ex)[i] <= this->x_lo_[i]+epsn && (*eg)[i] > 0.0) ){
530  (*ev)[i] = 0.0;
531  }
532  }
533  }
535  ROL::Ptr<const std::vector<Real> > ex =
536  dynamic_cast<const ROL::StdVector<Real>&>(x).getVector();
537  ROL::Ptr<const std::vector<Real> > eg =
538  dynamic_cast<const ROL::StdVector<Real>&>(g).getVector();
539  ROL::Ptr<std::vector<Real> > ev =
540  dynamic_cast<ROL::StdVector<Real>&>(v).getVector();
541  Real epsn = std::min(eps,this->min_diff_);
542  for ( int i = 0; i < this->dim_; i++ ) {
543  if ( ((*ex)[i] >= this->x_up_[i]-epsn && (*eg)[i] < 0.0) ) {
544  (*ev)[i] = 0.0;
545  }
546  }
547  }
549  ROL::Ptr<std::vector<Real> > us = ROL::makePtr<std::vector<Real>>(this->dim_,0.0);
550  us->assign(this->x_up_.begin(),this->x_up_.end());
551  ROL::Ptr<ROL::Vector<Real> > up = ROL::makePtr<ROL::StdVector<Real>>(us);
552  u.set(*up);
553  }
555  ROL::Ptr<std::vector<Real> > ls = ROL::makePtr<std::vector<Real>>(this->dim_,0.0);
556  ls->assign(this->x_lo_.begin(),this->x_lo_.end());
557  ROL::Ptr<ROL::Vector<Real> > lp = ROL::makePtr<ROL::StdVector<Real>>(ls);
558  l.set(*lp);
559  }
560 };
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:80
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 -binding set.
Defines a no-output stream class ROL::NullStream and a function makeStreamPtr which either wraps a re...
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)
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:209
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)
Real compute_norm(const std::vector< Real > &r)
void scale(std::vector< Real > &u, const Real alpha=0.0)