ROL
burgers-control/example_02.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"
52
53template<class Real>
55private:
56 int nx_;
57 Real dx_;
58 Real nu_;
59 Real u0_;
60 Real u1_;
61 Real f_;
62
63private:
64 Real compute_norm(const std::vector<Real> &r) {
65 return std::sqrt(this->dot(r,r));
66 }
67
68 Real dot(const std::vector<Real> &x, const std::vector<Real> &y) {
69 Real ip = 0.0;
70 Real c = (((int)x.size()==this->nx_) ? 4.0 : 2.0);
71 for (unsigned i=0; i<x.size(); i++) {
72 if ( i == 0 ) {
73 ip += this->dx_/6.0*(c*x[i] + x[i+1])*y[i];
74 }
75 else if ( i == x.size()-1 ) {
76 ip += this->dx_/6.0*(x[i-1] + c*x[i])*y[i];
77 }
78 else {
79 ip += this->dx_/6.0*(x[i-1] + 4.0*x[i] + x[i+1])*y[i];
80 }
81 }
82 return ip;
83 }
84
86
87 void update(std::vector<Real> &u, const std::vector<Real> &s, const Real alpha=1.0) {
88 for (unsigned i=0; i<u.size(); i++) {
89 u[i] += alpha*s[i];
90 }
91 }
92
93 void scale(std::vector<Real> &u, const Real alpha=0.0) {
94 for (unsigned i=0; i<u.size(); i++) {
95 u[i] *= alpha;
96 }
97 }
98
99 void compute_residual(std::vector<Real> &r, const std::vector<Real> &u,
100 const std::vector<Real> &z) {
101 r.clear();
102 r.resize(this->nx_,0.0);
103 for (int i=0; i<this->nx_; i++) {
104 // Contribution from stiffness term
105 if ( i==0 ) {
106 r[i] = this->nu_/this->dx_*(2.0*u[i]-u[i+1]);
107 }
108 else if (i==this->nx_-1) {
109 r[i] = this->nu_/this->dx_*(2.0*u[i]-u[i-1]);
110 }
111 else {
112 r[i] = this->nu_/this->dx_*(2.0*u[i]-u[i-1]-u[i+1]);
113 }
114 // Contribution from nonlinear term
115 if (i<this->nx_-1){
116 r[i] += u[i+1]*(u[i]+u[i+1])/6.0;
117 }
118 if (i>0) {
119 r[i] -= u[i-1]*(u[i-1]+u[i])/6.0;
120 }
121 // Contribution from control
122 r[i] -= this->dx_/6.0*(z[i]+4.0*z[i+1]+z[i+2]);
123 // Contribution from right-hand side
124 r[i] -= this->dx_*this->f_;
125 }
126 // Contribution from Dirichlet boundary terms
127 r[0] -= this->u0_*u[ 0]/6.0 + this->u0_*this->u0_/6.0 + this->nu_*this->u0_/this->dx_;
128 r[this->nx_-1] += this->u1_*u[this->nx_-1]/6.0 + this->u1_*this->u1_/6.0 - this->nu_*this->u1_/this->dx_;
129 }
130
131 void compute_pde_jacobian(std::vector<Real> &dl, std::vector<Real> &d, std::vector<Real> &du,
132 const std::vector<Real> &u) {
133 // Get Diagonal and Off-Diagonal Entries of linear PDE Jacobian
134 d.clear();
135 d.resize(this->nx_,this->nu_*2.0/this->dx_);
136 dl.clear();
137 dl.resize(this->nx_-1,-this->nu_/this->dx_);
138 du.clear();
139 du.resize(this->nx_-1,-this->nu_/this->dx_);
140 // Contribution from nonlinearity
141 for (int i=0; i<this->nx_; i++) {
142 if (i<this->nx_-1) {
143 dl[i] += (-2.0*u[i]-u[i+1])/6.0;
144 d[i] += u[i+1]/6.0;
145 }
146 if (i>0) {
147 d[i] += -u[i-1]/6.0;
148 du[i-1] += (u[i-1]+2.0*u[i])/6.0;
149 }
150 }
151 // Contribution from Dirichlet boundary conditions
152 d[0] -= this->u0_/6.0;
153 d[this->nx_-1] += this->u1_/6.0;
154 }
155
156 void linear_solve(std::vector<Real> &u, std::vector<Real> &dl, std::vector<Real> &d, std::vector<Real> &du,
157 const std::vector<Real> &r, const bool transpose = false) {
158 u.assign(r.begin(),r.end());
159 // Perform LDL factorization
160 Teuchos::LAPACK<int,Real> lp;
161 std::vector<Real> du2(this->nx_-2,0.0);
162 std::vector<int> ipiv(this->nx_,0);
163 int info;
164 int ldb = this->nx_;
165 int nhrs = 1;
166 lp.GTTRF(this->nx_,&dl[0],&d[0],&du[0],&du2[0],&ipiv[0],&info);
167 char trans = 'N';
168 if ( transpose ) {
169 trans = 'T';
170 }
171 lp.GTTRS(trans,this->nx_,nhrs,&dl[0],&d[0],&du[0],&du2[0],&ipiv[0],&u[0],ldb,&info);
172 }
173
174public:
175
176 Constraint_BurgersControl(int nx = 128, Real nu = 1.e-2, Real u0 = 1.0, Real u1 = 0.0, Real f = 0.0)
177 : nx_(nx), nu_(nu), u0_(u0), u1_(u1), f_(f) {
178 dx_ = 1.0/((Real)nx+1.0);
179 }
180
182 const ROL::Vector<Real> &z, Real &tol) {
183 ROL::Ptr<std::vector<Real> > cp =
184 dynamic_cast<ROL::StdVector<Real>&>(c).getVector();
185 ROL::Ptr<const std::vector<Real> > up =
186 dynamic_cast<const ROL::StdVector<Real>&>(u).getVector();
187 ROL::Ptr<const std::vector<Real> > zp =
188 dynamic_cast<const ROL::StdVector<Real>&>(z).getVector();
189 this->compute_residual(*cp,*up,*zp);
190 }
191
193 const ROL::Vector<Real> &z, Real &tol) {
194 ROL::Ptr<std::vector<Real> > jvp =
195 dynamic_cast<ROL::StdVector<Real>&>(jv).getVector();
196 ROL::Ptr<const std::vector<Real> > vp =
197 dynamic_cast<const ROL::StdVector<Real>&>(v).getVector();
198 ROL::Ptr<const std::vector<Real> > up =
199 dynamic_cast<const ROL::StdVector<Real>&>(u).getVector();
200 ROL::Ptr<const std::vector<Real> > zp =
201 dynamic_cast<const ROL::StdVector<Real>&>(z).getVector();
202 // Fill jvp
203 for (int i = 0; i < this->nx_; i++) {
204 (*jvp)[i] = this->nu_/this->dx_*2.0*(*vp)[i];
205 if ( i > 0 ) {
206 (*jvp)[i] += -this->nu_/this->dx_*(*vp)[i-1]
207 -(*up)[i-1]/6.0*(*vp)[i]
208 -((*up)[i]+2.0*(*up)[i-1])/6.0*(*vp)[i-1];
209 }
210 if ( i < this->nx_-1 ) {
211 (*jvp)[i] += -this->nu_/this->dx_*(*vp)[i+1]
212 +(*up)[i+1]/6.0*(*vp)[i]
213 +((*up)[i]+2.0*(*up)[i+1])/6.0*(*vp)[i+1];
214 }
215 }
216 (*jvp)[0] -= this->u0_/6.0*(*vp)[0];
217 (*jvp)[this->nx_-1] += this->u1_/6.0*(*vp)[this->nx_-1];
218 }
219
221 const ROL::Vector<Real> &z, Real &tol) {
222 ROL::Ptr<std::vector<Real> > jvp =
223 dynamic_cast<ROL::StdVector<Real>&>(jv).getVector();
224 ROL::Ptr<const std::vector<Real> > vp =
225 dynamic_cast<const ROL::StdVector<Real>&>(v).getVector();
226 ROL::Ptr<const std::vector<Real> > up =
227 dynamic_cast<const ROL::StdVector<Real>&>(u).getVector();
228 ROL::Ptr<const std::vector<Real> > zp =
229 dynamic_cast<const ROL::StdVector<Real>&>(z).getVector();
230 for (int i=0; i<this->nx_; i++) {
231 // Contribution from control
232 (*jvp)[i] = -this->dx_/6.0*((*vp)[i]+4.0*(*vp)[i+1]+(*vp)[i+2]);
233 }
234 }
235
237 const ROL::Vector<Real> &z, Real &tol) {
238 ROL::Ptr<std::vector<Real> > ijvp =
239 dynamic_cast<ROL::StdVector<Real>&>(ijv).getVector();
240 ROL::Ptr<const std::vector<Real> > vp =
241 dynamic_cast<const ROL::StdVector<Real>&>(v).getVector();
242 ROL::Ptr<const std::vector<Real> > up =
243 dynamic_cast<const ROL::StdVector<Real>&>(u).getVector();
244 ROL::Ptr<const std::vector<Real> > zp =
245 dynamic_cast<const ROL::StdVector<Real>&>(z).getVector();
246 // Get PDE Jacobian
247 std::vector<Real> d(this->nx_,0.0);
248 std::vector<Real> dl(this->nx_-1,0.0);
249 std::vector<Real> du(this->nx_-1,0.0);
250 this->compute_pde_jacobian(dl,d,du,*up);
251 // Solve solve state sensitivity system at current time step
252 this->linear_solve(*ijvp,dl,d,du,*vp);
253 }
254
256 const ROL::Vector<Real> &z, Real &tol) {
257 ROL::Ptr<std::vector<Real> > jvp =
258 dynamic_cast<ROL::StdVector<Real>&>(ajv).getVector();
259 ROL::Ptr<const std::vector<Real> > vp =
260 dynamic_cast<const ROL::StdVector<Real>&>(v).getVector();
261 ROL::Ptr<const std::vector<Real> > up =
262 dynamic_cast<const ROL::StdVector<Real>&>(u).getVector();
263 ROL::Ptr<const std::vector<Real> > zp =
264 dynamic_cast<const ROL::StdVector<Real>&>(z).getVector();
265 // Fill jvp
266 for (int i = 0; i < this->nx_; i++) {
267 (*jvp)[i] = this->nu_/this->dx_*2.0*(*vp)[i];
268 if ( i > 0 ) {
269 (*jvp)[i] += -this->nu_/this->dx_*(*vp)[i-1]
270 -(*up)[i-1]/6.0*(*vp)[i]
271 +((*up)[i-1]+2.0*(*up)[i])/6.0*(*vp)[i-1];
272 }
273 if ( i < this->nx_-1 ) {
274 (*jvp)[i] += -this->nu_/this->dx_*(*vp)[i+1]
275 +(*up)[i+1]/6.0*(*vp)[i]
276 -((*up)[i+1]+2.0*(*up)[i])/6.0*(*vp)[i+1];
277 }
278 }
279 (*jvp)[0] -= this->u0_/6.0*(*vp)[0];
280 (*jvp)[this->nx_-1] += this->u1_/6.0*(*vp)[this->nx_-1];
281 }
282
284 const ROL::Vector<Real> &z, Real &tol) {
285 ROL::Ptr<std::vector<Real> > jvp =
286 dynamic_cast<ROL::StdVector<Real>&>(jv).getVector();
287 ROL::Ptr<const std::vector<Real> > vp =
288 dynamic_cast<const ROL::StdVector<Real>&>(v).getVector();
289 ROL::Ptr<const std::vector<Real> > up =
290 dynamic_cast<const ROL::StdVector<Real>&>(u).getVector();
291 ROL::Ptr<const std::vector<Real> > zp =
292 dynamic_cast<const ROL::StdVector<Real>&>(z).getVector();
293 for (int i=0; i<this->nx_+2; i++) {
294 if ( i == 0 ) {
295 (*jvp)[i] = -this->dx_/6.0*(*vp)[i];
296 }
297 else if ( i == 1 ) {
298 (*jvp)[i] = -this->dx_/6.0*(4.0*(*vp)[i-1]+(*vp)[i]);
299 }
300 else if ( i == this->nx_ ) {
301 (*jvp)[i] = -this->dx_/6.0*(4.0*(*vp)[i-1]+(*vp)[i-2]);
302 }
303 else if ( i == this->nx_+1 ) {
304 (*jvp)[i] = -this->dx_/6.0*(*vp)[i-2];
305 }
306 else {
307 (*jvp)[i] = -this->dx_/6.0*((*vp)[i-2]+4.0*(*vp)[i-1]+(*vp)[i]);
308 }
309 }
310 }
311
313 const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
314 ROL::Ptr<std::vector<Real> > iajvp =
315 dynamic_cast<ROL::StdVector<Real>&>(iajv).getVector();
316 ROL::Ptr<const std::vector<Real> > vp =
317 dynamic_cast<const ROL::StdVector<Real>&>(v).getVector();
318 ROL::Ptr<const std::vector<Real> > up =
319 dynamic_cast<const ROL::StdVector<Real>&>(u).getVector();
320 // Get PDE Jacobian
321 std::vector<Real> d(this->nx_,0.0);
322 std::vector<Real> du(this->nx_-1,0.0);
323 std::vector<Real> dl(this->nx_-1,0.0);
324 this->compute_pde_jacobian(dl,d,du,*up);
325 // Solve solve adjoint system at current time step
326 this->linear_solve(*iajvp,dl,d,du,*vp,true);
327 }
328
330 const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
331 ROL::Ptr<std::vector<Real> > ahwvp =
332 dynamic_cast<ROL::StdVector<Real>&>(ahwv).getVector();
333 ROL::Ptr<const std::vector<Real> > wp =
334 dynamic_cast<const ROL::StdVector<Real>&>(w).getVector();
335 ROL::Ptr<const std::vector<Real> > vp =
336 dynamic_cast<const ROL::StdVector<Real>&>(v).getVector();
337 ROL::Ptr<const std::vector<Real> > up =
338 dynamic_cast<const ROL::StdVector<Real>&>(u).getVector();
339 ROL::Ptr<const std::vector<Real> > zp =
340 dynamic_cast<const ROL::StdVector<Real>&>(z).getVector();
341 for (int i=0; i<this->nx_; i++) {
342 // Contribution from nonlinear term
343 (*ahwvp)[i] = 0.0;
344 if (i<this->nx_-1){
345 (*ahwvp)[i] += ((*wp)[i]*(*vp)[i+1] - (*wp)[i+1]*(2.0*(*vp)[i]+(*vp)[i+1]))/6.0;
346 }
347 if (i>0) {
348 (*ahwvp)[i] += ((*wp)[i-1]*((*vp)[i-1]+2.0*(*vp)[i]) - (*wp)[i]*(*vp)[i-1])/6.0;
349 }
350 }
351 }
352
354 const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
355 ahwv.zero();
356 }
358 const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
359 ahwv.zero();
360 }
362 const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
363 ahwv.zero();
364 }
365
366// void solveAugmentedSystem(ROL::Vector<Real> &v1, ROL::Vector<Real> &v2, const ROL::Vector<Real> &b1,
367// const ROL::Vector<Real> &b2, const ROL::Vector<Real> &x, Real &tol) {}
368};
369
370template<class Real>
372private:
373 Real alpha_; // Penalty Parameter
374
375 int nx_; // Number of interior nodes
376 Real dx_; // Mesh spacing (i.e. 1/(nx+1))
377
378/***************************************************************/
379/********** BEGIN PRIVATE MEMBER FUNCTION DECLARATION **********/
380/***************************************************************/
381 Real evaluate_target(Real x) {
382 Real val = 0.0;
383 int example = 2;
384 switch (example) {
385 case 1: val = ((x<0.5) ? 1.0 : 0.0); break;
386 case 2: val = 1.0; break;
387 case 3: val = std::abs(std::sin(8.0*M_PI*x)); break;
388 case 4: val = std::exp(-0.5*(x-0.5)*(x-0.5)); break;
389 }
390 return val;
391 }
392
393 Real dot(const std::vector<Real> &x, const std::vector<Real> &y) {
394 Real ip = 0.0;
395 Real c = (((int)x.size()==this->nx_) ? 4.0 : 2.0);
396 for (unsigned i=0; i<x.size(); i++) {
397 if ( i == 0 ) {
398 ip += this->dx_/6.0*(c*x[i] + x[i+1])*y[i];
399 }
400 else if ( i == x.size()-1 ) {
401 ip += this->dx_/6.0*(x[i-1] + c*x[i])*y[i];
402 }
403 else {
404 ip += this->dx_/6.0*(x[i-1] + 4.0*x[i] + x[i+1])*y[i];
405 }
406 }
407 return ip;
408 }
409
410 void apply_mass(std::vector<Real> &Mu, const std::vector<Real> &u ) {
411 Mu.resize(u.size(),0.0);
412 Real c = (((int)u.size()==this->nx_) ? 4.0 : 2.0);
413 for (unsigned i=0; i<u.size(); i++) {
414 if ( i == 0 ) {
415 Mu[i] = this->dx_/6.0*(c*u[i] + u[i+1]);
416 }
417 else if ( i == u.size()-1 ) {
418 Mu[i] = this->dx_/6.0*(u[i-1] + c*u[i]);
419 }
420 else {
421 Mu[i] = this->dx_/6.0*(u[i-1] + 4.0*u[i] + u[i+1]);
422 }
423 }
424 }
425/*************************************************************/
426/********** END PRIVATE MEMBER FUNCTION DECLARATION **********/
427/*************************************************************/
428
429public:
430
431 Objective_BurgersControl(Real alpha = 1.e-4, int nx = 128) : alpha_(alpha), nx_(nx) {
432 dx_ = 1.0/((Real)nx+1.0);
433 }
434
435 Real value( const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
436 ROL::Ptr<const std::vector<Real> > up =
437 dynamic_cast<const ROL::StdVector<Real>&>(u).getVector();
438 ROL::Ptr<const std::vector<Real> > zp =
439 dynamic_cast<const ROL::StdVector<Real>&>(z).getVector();
440 // COMPUTE RESIDUAL
441 Real res1 = 0.0, res2 = 0.0, res3 = 0.0;
442 Real valu = 0.0, valz = this->dot(*zp,*zp);
443 for (int i=0; i<this->nx_; i++) {
444 if ( i == 0 ) {
445 res1 = (*up)[i]-evaluate_target((Real)(i+1)*this->dx_);
446 res2 = (*up)[i+1]-evaluate_target((Real)(i+2)*this->dx_);
447 valu += this->dx_/6.0*(4.0*res1 + res2)*res1;
448 }
449 else if ( i == this->nx_-1 ) {
450 res1 = (*up)[i-1]-evaluate_target((Real)i*this->dx_);
451 res2 = (*up)[i]-evaluate_target((Real)(i+1)*this->dx_);
452 valu += this->dx_/6.0*(res1 + 4.0*res2)*res2;
453 }
454 else {
455 res1 = (*up)[i-1]-evaluate_target((Real)i*this->dx_);
456 res2 = (*up)[i]-evaluate_target((Real)(i+1)*this->dx_);
457 res3 = (*up)[i+1]-evaluate_target((Real)(i+2)*this->dx_);
458 valu += this->dx_/6.0*(res1 + 4.0*res2 + res3)*res2;
459 }
460 }
461 return 0.5*(valu + this->alpha_*valz);
462 }
463
464 void gradient_1( ROL::Vector<Real> &g, const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
465 // Unwrap g
466 ROL::Ptr<std::vector<Real> > gup = ROL::constPtrCast<std::vector<Real> >(
467 (dynamic_cast<const ROL::StdVector<Real>&>(g)).getVector());
468 // Unwrap x
469 ROL::Ptr<const std::vector<Real> > up =
470 dynamic_cast<const ROL::StdVector<Real>&>(u).getVector();
471 ROL::Ptr<const std::vector<Real> > zp =
472 dynamic_cast<const ROL::StdVector<Real>&>(z).getVector();
473 // COMPUTE GRADIENT WRT U
474 std::vector<Real> diff(this->nx_,0.0);
475 for (int i=0; i<this->nx_; i++) {
476 diff[i] = ((*up)[i]-this->evaluate_target((Real)(i+1)*this->dx_));
477 }
478 this->apply_mass(*gup,diff);
479 }
480
481 void gradient_2( ROL::Vector<Real> &g, const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
482 // Unwrap g
483 ROL::Ptr<std::vector<Real> > gzp = ROL::constPtrCast<std::vector<Real> >(
484 (dynamic_cast<const ROL::StdVector<Real>&>(g)).getVector());
485 // Unwrap x
486 ROL::Ptr<const std::vector<Real> > up =
487 dynamic_cast<const ROL::StdVector<Real>&>(u).getVector();
488 ROL::Ptr<const std::vector<Real> > zp =
489 dynamic_cast<const ROL::StdVector<Real>&>(z).getVector();
490 // COMPUTE GRADIENT WRT Z
491 for (int i=0; i<this->nx_+2; i++) {
492 if (i==0) {
493 (*gzp)[i] = this->alpha_*this->dx_/6.0*(2.0*(*zp)[i]+(*zp)[i+1]);
494 }
495 else if (i==this->nx_+1) {
496 (*gzp)[i] = this->alpha_*this->dx_/6.0*(2.0*(*zp)[i]+(*zp)[i-1]);
497 }
498 else {
499 (*gzp)[i] = this->alpha_*this->dx_/6.0*((*zp)[i-1]+4.0*(*zp)[i]+(*zp)[i+1]);
500 }
501 }
502 }
503
505 const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
506 ROL::Ptr<std::vector<Real> > hvup =
507 dynamic_cast<ROL::StdVector<Real>&>(hv).getVector();
508 // Unwrap v
509 ROL::Ptr<const std::vector<Real> > vup =
510 dynamic_cast<const ROL::StdVector<Real>&>(v).getVector();
511 // COMPUTE GRADIENT WRT U
512 this->apply_mass(*hvup,*vup);
513 }
514
516 const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
517 hv.zero();
518 }
519
521 const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
522 hv.zero();
523 }
524
526 const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
527 ROL::Ptr<std::vector<Real> > hvzp =
528 dynamic_cast<ROL::StdVector<Real>&>(hv).getVector();
529 // Unwrap v
530 ROL::Ptr<const std::vector<Real> > vzp =
531 dynamic_cast<const ROL::StdVector<Real>&>(v).getVector();
532 // COMPUTE GRADIENT WRT Z
533 for (int i=0; i<this->nx_+2; i++) {
534 if (i==0) {
535 (*hvzp)[i] = this->alpha_*this->dx_/6.0*(2.0*(*vzp)[i]+(*vzp)[i+1]);
536 }
537 else if (i==this->nx_+1) {
538 (*hvzp)[i] = this->alpha_*this->dx_/6.0*(2.0*(*vzp)[i]+(*vzp)[i-1]);
539 }
540 else {
541 (*hvzp)[i] = this->alpha_*this->dx_/6.0*((*vzp)[i-1]+4.0*(*vzp)[i]+(*vzp)[i+1]);
542 }
543 }
544 }
545};
void scale(std::vector< Real > &u, const Real alpha=0.0)
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 ...
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 compute_residual(std::vector< Real > &r, const std::vector< Real > &u, const std::vector< Real > &z)
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...
Constraint_BurgersControl(int nx=128, Real nu=1.e-2, Real u0=1.0, Real u1=0.0, Real f=0.0)
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 compute_pde_jacobian(std::vector< Real > &dl, std::vector< Real > &d, std::vector< Real > &du, const std::vector< Real > &u)
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 .
Real compute_norm(const std::vector< Real > &r)
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...
void update(std::vector< Real > &u, const std::vector< Real > &s, const Real alpha=1.0)
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 .
void value(ROL::Vector< Real > &c, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Evaluate the constraint operator at .
Real dot(const std::vector< Real > &x, const std::vector< Real > &y)
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...
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)
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 .
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 .
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 hessVec_21(ROL::Vector< Real > &hv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
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.
Objective_BurgersControl(Real alpha=1.e-4, int nx=128)
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 hessVec_22(ROL::Vector< Real > &hv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
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.
void hessVec_12(ROL::Vector< Real > &hv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Real value(const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Compute value.
Defines the constraint operator interface for simulation-based optimization.
Provides the interface to evaluate simulation-based 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 zero()
Set to zero vector.