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