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