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 
33 template<class Real>
34 class Constraint_BurgersControl : public ROL::Constraint_SimOpt<Real> {
35 private:
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 
156 public:
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 
363 template<class Real>
364 class Objective_BurgersControl : public ROL::Objective_SimOpt<Real> {
365 private:
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 
422 public:
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::Objective_SimOpt
Provides the interface to evaluate simulation-based objective functions.
Definition: ROL_Objective_SimOpt.hpp:25
Objective_BurgersControl::dx_
Real dx_
Definition: burgers-control/example_01.hpp:25
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
ROL_PrimalDualRisk.hpp
ROL_Reduced_Objective_SimOpt.hpp
Objective_BurgersControl::evaluate_target
Real evaluate_target(Real x)
Definition: example_10.hpp:374
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
ROL_OptimizationSolver.hpp
Constraint_BurgersControl
Definition: test_04.hpp:670
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
ROL_Types.hpp
Contains definitions of custom data types in ROL.
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
ROL::Constraint::getParameter
const std::vector< Real > getParameter(void) const
Definition: ROL_Constraint.hpp:369
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
Objective_BurgersControl::nx_
int nx_
Definition: burgers-control/example_01.hpp:24
ROL::Vector::zero
virtual void zero()
Set to zero vector.
Definition: ROL_Vector.hpp:133
ROL::Vector
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:46
ROL_Stream.hpp
Defines a no-output stream class ROL::NullStream and a function makeStreamPtr which either wraps a re...
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::value
Real value(const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Compute value.
Definition: example_10.hpp:430
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
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::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::Constraint_BurgersControl
Constraint_BurgersControl(int nx=128)
Definition: example_10.hpp:158
Objective_BurgersControl::dot
Real dot(const std::vector< Real > &x, const std::vector< Real > &y)
Definition: example_10.hpp:386
ROL_StdVector.hpp
Objective_BurgersControl::Objective_BurgersControl
Objective_BurgersControl(Real alpha=1.e-4, int nx=128)
Definition: example_10.hpp:424
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
ROL_Constraint_SimOpt.hpp
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::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
Definition: burgers-control/example_01.hpp:20
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
Constraint_BurgersControl::dot
Real dot(const std::vector< Real > &x, const std::vector< Real > &y)
Definition: burgers-control/example_02.hpp:34
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::compute_residual
void compute_residual(std::vector< Real > &r, const std::vector< Real > &u, const std::vector< Real > &z)
Definition: example_10.hpp:74
Objective_BurgersControl::alpha_
Real alpha_
Definition: burgers-control/example_01.hpp:22
Constraint_BurgersControl::scale
void scale(std::vector< Real > &u, const Real alpha=0.0)
Definition: example_10.hpp:68
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
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
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
Objective_BurgersControl::apply_mass
void apply_mass(std::vector< Real > &Mu, const std::vector< Real > &u)
Definition: example_10.hpp:403
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::compute_norm
Real compute_norm(const std::vector< Real > &r)
Definition: example_10.hpp:39
Constraint_BurgersControl::nx_
int nx_
Definition: burgers-control/example_02.hpp:22
Constraint_BurgersControl::update
void update(std::vector< Real > &u, const std::vector< Real > &s, const Real alpha=1.0)
Definition: example_10.hpp:62
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
ROL::Constraint_SimOpt
Defines the constraint operator interface for simulation-based optimization.
Definition: ROL_Constraint_SimOpt.hpp:21
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::dx_
Real dx_
Definition: burgers-control/example_02.hpp:23
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
ROL_MonteCarloGenerator.hpp
ROL_Objective_SimOpt.hpp
Generated by   doxygen 1.8.5