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