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