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