ROL
example_05.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ************************************************************************
3 //
4 // Rapid Optimization Library (ROL) Package
5 // Copyright (2014) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact lead developers:
38 // Drew Kouri (dpkouri@sandia.gov) and
39 // Denis Ridzal (dridzal@sandia.gov)
40 //
41 // ************************************************************************
42 // @HEADER
43 
44 #include "ROL_Stream.hpp"
45 #include "Teuchos_LAPACK.hpp"
46 #include "Teuchos_GlobalMPISession.hpp"
47 #include "Teuchos_Comm.hpp"
48 #include "Teuchos_DefaultComm.hpp"
49 #include "Teuchos_CommHelpers.hpp"
50 
51 #include "ROL_ParameterList.hpp"
52 // ROL_Types contains predefined constants and objects
53 #include "ROL_Types.hpp"
54 // ROL vectors
55 #include "ROL_StdVector.hpp"
56 // ROL objective functions and constraints
57 #include "ROL_Objective_SimOpt.hpp"
58 #include "ROL_Constraint_SimOpt.hpp"
59 #include "ROL_Reduced_Objective_SimOpt.hpp"
60 #include "ROL_OptimizationSolver.hpp"
61 // ROL sample generators
62 #include "ROL_MonteCarloGenerator.hpp"
63 #include "ROL_StdTeuchosBatchManager.hpp"
64 
65 template<class Real>
66 class Constraint_BurgersControl : public ROL::Constraint_SimOpt<Real> {
67 private:
68  int nx_;
69  Real dx_;
70 
71  Real compute_norm(const std::vector<Real> &r) {
72  return std::sqrt(dot(r,r));
73  }
74 
75  Real dot(const std::vector<Real> &x, const std::vector<Real> &y) {
76  Real ip = 0.0;
77  Real c = (((int)x.size()==nx_) ? 4.0 : 2.0);
78  for (unsigned i=0; i<x.size(); i++) {
79  if ( i == 0 ) {
80  ip += dx_/6.0*(c*x[i] + x[i+1])*y[i];
81  }
82  else if ( i == x.size()-1 ) {
83  ip += dx_/6.0*(x[i-1] + c*x[i])*y[i];
84  }
85  else {
86  ip += dx_/6.0*(x[i-1] + 4.0*x[i] + x[i+1])*y[i];
87  }
88  }
89  return ip;
90  }
91 
92  using ROL::Constraint_SimOpt<Real>::update;
93 
94  void update(std::vector<Real> &u, const std::vector<Real> &s, const Real alpha=1.0) {
95  for (unsigned i=0; i<u.size(); i++) {
96  u[i] += alpha*s[i];
97  }
98  }
99 
100  void scale(std::vector<Real> &u, const Real alpha=0.0) {
101  for (unsigned i=0; i<u.size(); i++) {
102  u[i] *= alpha;
103  }
104  }
105 
106  void compute_residual(std::vector<Real> &r, const std::vector<Real> &u,
107  const std::vector<Real> &z) {
108  r.clear(); r.resize(nx_,0.0);
109  const std::vector<Real> param =
110  ROL::Constraint_SimOpt<Real>::getParameter();
111  Real nu = std::pow(10.0,param[0]-2.0);
112  Real f = param[1]/100.0;
113  Real u0 = 1.0+param[2]/1000.0;
114  Real u1 = param[3]/1000.0;
115  for (int i=0; i<nx_; i++) {
116  // Contribution from stiffness term
117  if ( i==0 ) {
118  r[i] = nu/dx_*(2.0*u[i]-u[i+1]);
119  }
120  else if (i==nx_-1) {
121  r[i] = nu/dx_*(2.0*u[i]-u[i-1]);
122  }
123  else {
124  r[i] = nu/dx_*(2.0*u[i]-u[i-1]-u[i+1]);
125  }
126  // Contribution from nonlinear term
127  if (i<nx_-1){
128  r[i] += u[i+1]*(u[i]+u[i+1])/6.0;
129  }
130  if (i>0) {
131  r[i] -= u[i-1]*(u[i-1]+u[i])/6.0;
132  }
133  // Contribution from control
134  r[i] -= dx_/6.0*(z[i]+4.0*z[i+1]+z[i+2]);
135  // Contribution from right-hand side
136  r[i] -= dx_*f;
137  }
138  // Contribution from Dirichlet boundary terms
139  r[ 0] -= u0*u[ 0]/6.0 + u0*u0/6.0 + nu*u0/dx_;
140  r[nx_-1] += u1*u[nx_-1]/6.0 + u1*u1/6.0 - nu*u1/dx_;
141  }
142 
143  void compute_pde_jacobian(std::vector<Real> &dl, std::vector<Real> &d, std::vector<Real> &du,
144  const std::vector<Real> &u) {
145  const std::vector<Real> param =
146  ROL::Constraint_SimOpt<Real>::getParameter();
147  Real nu = std::pow(10.0,param[0]-2.0);
148  Real u0 = 1.0+param[2]/1000.0;
149  Real u1 = param[3]/1000.0;
150  // Get Diagonal and Off-Diagonal Entries of linear PDE Jacobian
151  d.clear(); d.resize(nx_,nu*2.0/dx_);
152  dl.clear(); dl.resize(nx_-1,-nu/dx_);
153  du.clear(); du.resize(nx_-1,-nu/dx_);
154  // Contribution from nonlinearity
155  for (int i=0; i<nx_; i++) {
156  if (i<nx_-1) {
157  dl[i] += (-2.0*u[i]-u[i+1])/6.0;
158  d[i] += u[i+1]/6.0;
159  }
160  if (i>0) {
161  d[i] += -u[i-1]/6.0;
162  du[i-1] += (u[i-1]+2.0*u[i])/6.0;
163  }
164  }
165  // Contribution from Dirichlet boundary conditions
166  d[ 0] -= u0/6.0;
167  d[nx_-1] += u1/6.0;
168  }
169 
170  void linear_solve(std::vector<Real> &u, std::vector<Real> &dl, std::vector<Real> &d, std::vector<Real> &du,
171  const std::vector<Real> &r, const bool transpose = false) {
172  u.assign(r.begin(),r.end());
173  // Perform LDL factorization
174  Teuchos::LAPACK<int,Real> lp;
175  std::vector<Real> du2(nx_-2,0.0);
176  std::vector<int> ipiv(nx_,0);
177  int info;
178  int ldb = nx_;
179  int nhrs = 1;
180  lp.GTTRF(nx_,&dl[0],&d[0],&du[0],&du2[0],&ipiv[0],&info);
181  char trans = 'N';
182  if ( transpose ) {
183  trans = 'T';
184  }
185  lp.GTTRS(trans,nx_,nhrs,&dl[0],&d[0],&du[0],&du2[0],&ipiv[0],&u[0],ldb,&info);
186  }
187 
188 public:
189 
190  Constraint_BurgersControl(int nx = 128) : nx_(nx), dx_(1.0/((Real)nx+1.0)) {}
191 
192  void value(ROL::Vector<Real> &c, const ROL::Vector<Real> &u,
193  const ROL::Vector<Real> &z, Real &tol) {
194  ROL::Ptr<std::vector<Real> > cp =
195  dynamic_cast<ROL::StdVector<Real>&>(c).getVector();
196  ROL::Ptr<const std::vector<Real> > up =
197  dynamic_cast<const ROL::StdVector<Real>&>(u).getVector();
198  ROL::Ptr<const std::vector<Real> > zp =
199  dynamic_cast<const ROL::StdVector<Real>&>(z).getVector();
200  compute_residual(*cp,*up,*zp);
201  }
202 
203  void solve(ROL::Vector<Real> &c, ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
204  ROL::Ptr<std::vector<Real> > up =
205  ROL::constPtrCast<std::vector<Real> >((dynamic_cast<ROL::StdVector<Real>&>(u)).getVector());
206  up->assign(up->size(),static_cast<Real>(1));
207  ROL::Constraint_SimOpt<Real>::solve(c,u,z,tol);
208  }
209 
210  void applyJacobian_1(ROL::Vector<Real> &jv, const ROL::Vector<Real> &v, const ROL::Vector<Real> &u,
211  const ROL::Vector<Real> &z, Real &tol) {
212  ROL::Ptr<std::vector<Real> > jvp =
213  dynamic_cast<ROL::StdVector<Real>&>(jv).getVector();
214  ROL::Ptr<const std::vector<Real> > vp =
215  dynamic_cast<const ROL::StdVector<Real>&>(v).getVector();
216  ROL::Ptr<const std::vector<Real> > up =
217  dynamic_cast<const ROL::StdVector<Real>&>(u).getVector();
218  ROL::Ptr<const std::vector<Real> > zp =
219  dynamic_cast<const ROL::StdVector<Real>&>(z).getVector();
220  const std::vector<Real> param =
221  ROL::Constraint_SimOpt<Real>::getParameter();
222  Real nu = std::pow(10.0,param[0]-2.0);
223  Real u0 = 1.0+param[2]/1000.0;
224  Real u1 = param[3]/1000.0;
225  // Fill jvp
226  for (int i = 0; i < nx_; i++) {
227  (*jvp)[i] = nu/dx_*2.0*(*vp)[i];
228  if ( i > 0 ) {
229  (*jvp)[i] += -nu/dx_*(*vp)[i-1]
230  -(*up)[i-1]/6.0*(*vp)[i]
231  -((*up)[i]+2.0*(*up)[i-1])/6.0*(*vp)[i-1];
232  }
233  if ( i < nx_-1 ) {
234  (*jvp)[i] += -nu/dx_*(*vp)[i+1]
235  +(*up)[i+1]/6.0*(*vp)[i]
236  +((*up)[i]+2.0*(*up)[i+1])/6.0*(*vp)[i+1];
237  }
238  }
239  (*jvp)[ 0] -= u0/6.0*(*vp)[0];
240  (*jvp)[nx_-1] += u1/6.0*(*vp)[nx_-1];
241  }
242 
243  void applyJacobian_2(ROL::Vector<Real> &jv, const ROL::Vector<Real> &v, const ROL::Vector<Real> &u,
244  const ROL::Vector<Real> &z, Real &tol) {
245  ROL::Ptr<std::vector<Real> > jvp =
246  dynamic_cast<ROL::StdVector<Real>&>(jv).getVector();
247  ROL::Ptr<const std::vector<Real>> vp =
248  dynamic_cast<const ROL::StdVector<Real>&>(v).getVector();
249  ROL::Ptr<const std::vector<Real>> up =
250  dynamic_cast<const ROL::StdVector<Real>&>(u).getVector();
251  ROL::Ptr<const std::vector<Real>> zp =
252  dynamic_cast<const ROL::StdVector<Real>&>(z).getVector();
253  for (int i=0; i<nx_; i++) {
254  // Contribution from control
255  (*jvp)[i] = -dx_/6.0*((*vp)[i]+4.0*(*vp)[i+1]+(*vp)[i+2]);
256  }
257  }
258 
259  void applyInverseJacobian_1(ROL::Vector<Real> &ijv, const ROL::Vector<Real> &v, const ROL::Vector<Real> &u,
260  const ROL::Vector<Real> &z, Real &tol) {
261  ROL::Ptr<std::vector<Real> > ijvp =
262  dynamic_cast<ROL::StdVector<Real>&>(ijv).getVector();
263  ROL::Ptr<const std::vector<Real> > vp =
264  dynamic_cast<const ROL::StdVector<Real>&>(v).getVector();
265  ROL::Ptr<const std::vector<Real> > up =
266  dynamic_cast<const ROL::StdVector<Real>&>(u).getVector();
267  ROL::Ptr<const std::vector<Real> > zp =
268  dynamic_cast<const ROL::StdVector<Real>&>(z).getVector();
269  // Get PDE Jacobian
270  std::vector<Real> d(nx_,0.0);
271  std::vector<Real> dl(nx_-1,0.0);
272  std::vector<Real> du(nx_-1,0.0);
273  compute_pde_jacobian(dl,d,du,*up);
274  // Solve solve state sensitivity system at current time step
275  linear_solve(*ijvp,dl,d,du,*vp);
276  }
277 
278  void applyAdjointJacobian_1(ROL::Vector<Real> &ajv, const ROL::Vector<Real> &v, const ROL::Vector<Real> &u,
279  const ROL::Vector<Real> &z, Real &tol) {
280  ROL::Ptr<std::vector<Real> > jvp =
281  dynamic_cast<ROL::StdVector<Real>&>(ajv).getVector();
282  ROL::Ptr<const std::vector<Real> > vp =
283  dynamic_cast<const ROL::StdVector<Real>&>(v).getVector();
284  ROL::Ptr<const std::vector<Real> > up =
285  dynamic_cast<const ROL::StdVector<Real>&>(u).getVector();
286  ROL::Ptr<const std::vector<Real> > zp =
287  dynamic_cast<const ROL::StdVector<Real>&>(z).getVector();
288  const std::vector<Real> param =
289  ROL::Constraint_SimOpt<Real>::getParameter();
290  Real nu = std::pow(10.0,param[0]-2.0);
291  Real u0 = 1.0+param[2]/1000.0;
292  Real u1 = param[3]/1000.0;
293  // Fill jvp
294  for (int i = 0; i < nx_; i++) {
295  (*jvp)[i] = nu/dx_*2.0*(*vp)[i];
296  if ( i > 0 ) {
297  (*jvp)[i] += -nu/dx_*(*vp)[i-1]
298  -(*up)[i-1]/6.0*(*vp)[i]
299  +((*up)[i-1]+2.0*(*up)[i])/6.0*(*vp)[i-1];
300  }
301  if ( i < nx_-1 ) {
302  (*jvp)[i] += -nu/dx_*(*vp)[i+1]
303  +(*up)[i+1]/6.0*(*vp)[i]
304  -((*up)[i+1]+2.0*(*up)[i])/6.0*(*vp)[i+1];
305  }
306  }
307  (*jvp)[ 0] -= u0/6.0*(*vp)[0];
308  (*jvp)[nx_-1] += u1/6.0*(*vp)[nx_-1];
309  }
310 
311  void applyAdjointJacobian_2(ROL::Vector<Real> &jv, const ROL::Vector<Real> &v, const ROL::Vector<Real> &u,
312  const ROL::Vector<Real> &z, Real &tol) {
313  ROL::Ptr<std::vector<Real> > jvp =
314  dynamic_cast<ROL::StdVector<Real>&>(jv).getVector();
315  ROL::Ptr<const std::vector<Real> > vp =
316  dynamic_cast<const ROL::StdVector<Real>&>(v).getVector();
317  ROL::Ptr<const std::vector<Real> > up =
318  dynamic_cast<const ROL::StdVector<Real>&>(u).getVector();
319  ROL::Ptr<const std::vector<Real> > zp =
320  dynamic_cast<const ROL::StdVector<Real>&>(z).getVector();
321  for (int i=0; i<nx_+2; i++) {
322  if ( i == 0 ) {
323  (*jvp)[i] = -dx_/6.0*(*vp)[i];
324  }
325  else if ( i == 1 ) {
326  (*jvp)[i] = -dx_/6.0*(4.0*(*vp)[i-1]+(*vp)[i]);
327  }
328  else if ( i == nx_ ) {
329  (*jvp)[i] = -dx_/6.0*(4.0*(*vp)[i-1]+(*vp)[i-2]);
330  }
331  else if ( i == nx_+1 ) {
332  (*jvp)[i] = -dx_/6.0*(*vp)[i-2];
333  }
334  else {
335  (*jvp)[i] = -dx_/6.0*((*vp)[i-2]+4.0*(*vp)[i-1]+(*vp)[i]);
336  }
337  }
338  }
339 
340  void applyInverseAdjointJacobian_1(ROL::Vector<Real> &iajv, const ROL::Vector<Real> &v,
341  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
342  ROL::Ptr<std::vector<Real> > iajvp =
343  dynamic_cast<ROL::StdVector<Real>&>(iajv).getVector();
344  ROL::Ptr<const std::vector<Real> > vp =
345  dynamic_cast<const ROL::StdVector<Real>&>(v).getVector();
346  ROL::Ptr<const std::vector<Real>> up =
347  dynamic_cast<const ROL::StdVector<Real>&>(u).getVector();
348  // Get PDE Jacobian
349  std::vector<Real> d(nx_,0.0);
350  std::vector<Real> du(nx_-1,0.0);
351  std::vector<Real> dl(nx_-1,0.0);
352  compute_pde_jacobian(dl,d,du,*up);
353  // Solve solve adjoint system at current time step
354  linear_solve(*iajvp,dl,d,du,*vp,true);
355  }
356 
357  void applyAdjointHessian_11(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  ROL::Ptr<std::vector<Real> > ahwvp =
360  dynamic_cast<ROL::StdVector<Real>&>(ahwv).getVector();
361  ROL::Ptr<const std::vector<Real> > wp =
362  dynamic_cast<const ROL::StdVector<Real>&>(w).getVector();
363  ROL::Ptr<const std::vector<Real> > vp =
364  dynamic_cast<const ROL::StdVector<Real>&>(v).getVector();
365  ROL::Ptr<const std::vector<Real> > up =
366  dynamic_cast<const ROL::StdVector<Real>&>(u).getVector();
367  ROL::Ptr<const std::vector<Real> > zp =
368  dynamic_cast<const ROL::StdVector<Real>&>(z).getVector();
369  for (int i=0; i<nx_; i++) {
370  // Contribution from nonlinear term
371  (*ahwvp)[i] = 0.0;
372  if (i<nx_-1){
373  (*ahwvp)[i] += ((*wp)[i]*(*vp)[i+1] - (*wp)[i+1]*(2.0*(*vp)[i]+(*vp)[i+1]))/6.0;
374  }
375  if (i>0) {
376  (*ahwvp)[i] += ((*wp)[i-1]*((*vp)[i-1]+2.0*(*vp)[i]) - (*wp)[i]*(*vp)[i-1])/6.0;
377  }
378  }
379  }
380 
381  void applyAdjointHessian_12(ROL::Vector<Real> &ahwv, const ROL::Vector<Real> &w, const ROL::Vector<Real> &v,
382  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
383  ahwv.zero();
384  }
385  void applyAdjointHessian_21(ROL::Vector<Real> &ahwv, const ROL::Vector<Real> &w, const ROL::Vector<Real> &v,
386  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
387  ahwv.zero();
388  }
389  void applyAdjointHessian_22(ROL::Vector<Real> &ahwv, const ROL::Vector<Real> &w, const ROL::Vector<Real> &v,
390  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
391  ahwv.zero();
392  }
393 };
394 
395 template<class Real>
396 class Objective_BurgersControl : public ROL::Objective_SimOpt<Real> {
397 private:
398  Real alpha_; // Penalty Parameter
399 
400  int nx_; // Number of interior nodes
401  Real dx_; // Mesh spacing (i.e. 1/(nx+1))
402 
403 /***************************************************************/
404 /********** BEGIN PRIVATE MEMBER FUNCTION DECLARATION **********/
405 /***************************************************************/
406  Real evaluate_target(Real x) {
407  Real val = 0.0;
408  int example = 2;
409  switch (example) {
410  case 1: val = ((x<0.5) ? 1.0 : 0.0); break;
411  case 2: val = 1.0; break;
412  case 3: val = std::abs(std::sin(8.0*M_PI*x)); break;
413  case 4: val = std::exp(-0.5*(x-0.5)*(x-0.5)); break;
414  }
415  return val;
416  }
417 
418  Real dot(const std::vector<Real> &x, const std::vector<Real> &y) {
419  Real ip = 0.0;
420  Real c = (((int)x.size()==nx_) ? 4.0 : 2.0);
421  for (unsigned i=0; i<x.size(); i++) {
422  if ( i == 0 ) {
423  ip += dx_/6.0*(c*x[i] + x[i+1])*y[i];
424  }
425  else if ( i == x.size()-1 ) {
426  ip += dx_/6.0*(x[i-1] + c*x[i])*y[i];
427  }
428  else {
429  ip += dx_/6.0*(x[i-1] + 4.0*x[i] + x[i+1])*y[i];
430  }
431  }
432  return ip;
433  }
434 
435  void apply_mass(std::vector<Real> &Mu, const std::vector<Real> &u ) {
436  Mu.resize(u.size(),0.0);
437  Real c = (((int)u.size()==nx_) ? 4.0 : 2.0);
438  for (unsigned i=0; i<u.size(); i++) {
439  if ( i == 0 ) {
440  Mu[i] = dx_/6.0*(c*u[i] + u[i+1]);
441  }
442  else if ( i == u.size()-1 ) {
443  Mu[i] = dx_/6.0*(u[i-1] + c*u[i]);
444  }
445  else {
446  Mu[i] = dx_/6.0*(u[i-1] + 4.0*u[i] + u[i+1]);
447  }
448  }
449  }
450 /*************************************************************/
451 /********** END PRIVATE MEMBER FUNCTION DECLARATION **********/
452 /*************************************************************/
453 
454 public:
455 
456  Objective_BurgersControl(Real alpha = 1.e-4, int nx = 128) : alpha_(alpha), nx_(nx) {
457  dx_ = 1.0/((Real)nx+1.0);
458  }
459 
460  using ROL::Objective_SimOpt<Real>::value;
461 
462  Real value( const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
463  ROL::Ptr<const std::vector<Real> > up =
464  dynamic_cast<const ROL::StdVector<Real>&>(u).getVector();
465  ROL::Ptr<const std::vector<Real> > zp =
466  dynamic_cast<const ROL::StdVector<Real>&>(z).getVector();
467  // COMPUTE RESIDUAL
468  Real res1 = 0.0, res2 = 0.0, res3 = 0.0;
469  Real valu = 0.0, valz = dot(*zp,*zp);
470  for (int i=0; i<nx_; i++) {
471  if ( i == 0 ) {
472  res1 = (*up)[i]-evaluate_target((Real)(i+1)*dx_);
473  res2 = (*up)[i+1]-evaluate_target((Real)(i+2)*dx_);
474  valu += dx_/6.0*(4.0*res1 + res2)*res1;
475  }
476  else if ( i == nx_-1 ) {
477  res1 = (*up)[i-1]-evaluate_target((Real)i*dx_);
478  res2 = (*up)[i]-evaluate_target((Real)(i+1)*dx_);
479  valu += dx_/6.0*(res1 + 4.0*res2)*res2;
480  }
481  else {
482  res1 = (*up)[i-1]-evaluate_target((Real)i*dx_);
483  res2 = (*up)[i]-evaluate_target((Real)(i+1)*dx_);
484  res3 = (*up)[i+1]-evaluate_target((Real)(i+2)*dx_);
485  valu += dx_/6.0*(res1 + 4.0*res2 + res3)*res2;
486  }
487  }
488  return 0.5*(valu + alpha_*valz);
489  }
490 
491  void gradient_1( ROL::Vector<Real> &g, const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
492  // Unwrap g
493  ROL::Ptr<std::vector<Real> > gup =
494  dynamic_cast<ROL::StdVector<Real>&>(g).getVector();
495  // Unwrap x
496  ROL::Ptr<const std::vector<Real> > up =
497  dynamic_cast<const ROL::StdVector<Real>&>(u).getVector();
498  ROL::Ptr<const std::vector<Real> > zp =
499  dynamic_cast<const ROL::StdVector<Real>&>(z).getVector();
500  // COMPUTE GRADIENT WRT U
501  std::vector<Real> diff(nx_,0.0);
502  for (int i=0; i<nx_; i++) {
503  diff[i] = ((*up)[i]-evaluate_target((Real)(i+1)*dx_));
504  }
505  apply_mass(*gup,diff);
506  }
507 
508  void gradient_2( ROL::Vector<Real> &g, const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
509  // Unwrap g
510  ROL::Ptr<std::vector<Real> > gzp =
511  dynamic_cast<ROL::StdVector<Real>&>(g).getVector();
512  // Unwrap x
513  ROL::Ptr<const std::vector<Real> > up =
514  dynamic_cast<const ROL::StdVector<Real>&>(u).getVector();
515  ROL::Ptr<const std::vector<Real> > zp =
516  dynamic_cast<const ROL::StdVector<Real>&>(z).getVector();
517  // COMPUTE GRADIENT WRT Z
518  for (int i=0; i<nx_+2; i++) {
519  if (i==0) {
520  (*gzp)[i] = alpha_*dx_/6.0*(2.0*(*zp)[i]+(*zp)[i+1]);
521  }
522  else if (i==nx_+1) {
523  (*gzp)[i] = alpha_*dx_/6.0*(2.0*(*zp)[i]+(*zp)[i-1]);
524  }
525  else {
526  (*gzp)[i] = alpha_*dx_/6.0*((*zp)[i-1]+4.0*(*zp)[i]+(*zp)[i+1]);
527  }
528  }
529  }
530 
531  void hessVec_11( ROL::Vector<Real> &hv, const ROL::Vector<Real> &v,
532  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
533  ROL::Ptr<std::vector<Real> > hvup =
534  dynamic_cast<ROL::StdVector<Real>&>(hv).getVector();
535  // Unwrap v
536  ROL::Ptr<const std::vector<Real> > vup =
537  dynamic_cast<const ROL::StdVector<Real>&>(v).getVector();
538  // COMPUTE GRADIENT WRT U
539  apply_mass(*hvup,*vup);
540  }
541 
542  void hessVec_12( ROL::Vector<Real> &hv, const ROL::Vector<Real> &v,
543  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
544  hv.zero();
545  }
546 
547  void hessVec_21( ROL::Vector<Real> &hv, const ROL::Vector<Real> &v,
548  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
549  hv.zero();
550  }
551 
552  void hessVec_22( ROL::Vector<Real> &hv, const ROL::Vector<Real> &v,
553  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
554  ROL::Ptr<std::vector<Real> > hvzp =
555  dynamic_cast<ROL::StdVector<Real>&>(hv).getVector();
556  // Unwrap v
557  ROL::Ptr<const std::vector<Real> > vzp =
558  dynamic_cast<const ROL::StdVector<Real>&>(v).getVector();
559  // COMPUTE GRADIENT WRT Z
560  for (int i=0; i<nx_+2; i++) {
561  if (i==0) {
562  (*hvzp)[i] = alpha_*dx_/6.0*(2.0*(*vzp)[i]+(*vzp)[i+1]);
563  }
564  else if (i==nx_+1) {
565  (*hvzp)[i] = alpha_*dx_/6.0*(2.0*(*vzp)[i]+(*vzp)[i-1]);
566  }
567  else {
568  (*hvzp)[i] = alpha_*dx_/6.0*((*vzp)[i-1]+4.0*(*vzp)[i]+(*vzp)[i+1]);
569  }
570  }
571  }
572 };
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:72
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:278
ROL_Reduced_Objective_SimOpt.hpp
Objective_BurgersControl::evaluate_target
Real evaluate_target(Real x)
Definition: example_05.hpp:406
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:381
ROL_OptimizationSolver.hpp
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:203
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:357
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:340
Objective_BurgersControl::nx_
int nx_
Definition: burgers-control/example_01.hpp:71
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:552
Objective_BurgersControl::value
Real value(const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Compute value.
Definition: example_05.hpp:462
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:491
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:385
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:143
Constraint_BurgersControl::Constraint_BurgersControl
Constraint_BurgersControl(int nx=128)
Definition: example_05.hpp:190
Objective_BurgersControl::dot
Real dot(const std::vector< Real > &x, const std::vector< Real > &y)
Definition: example_05.hpp:418
ROL_StdVector.hpp
Objective_BurgersControl::Objective_BurgersControl
Objective_BurgersControl(Real alpha=1.e-4, int nx=128)
Definition: example_05.hpp:456
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:210
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:547
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:542
Objective_BurgersControl
Definition: burgers-control/example_01.hpp:67
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:508
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:259
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:106
Objective_BurgersControl::alpha_
Real alpha_
Definition: burgers-control/example_01.hpp:69
Constraint_BurgersControl::scale
void scale(std::vector< Real > &u, const Real alpha=0.0)
Definition: example_05.hpp:100
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:531
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:243
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:207
Objective_BurgersControl::apply_mass
void apply_mass(std::vector< Real > &Mu, const std::vector< Real > &u)
Definition: example_05.hpp:435
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:389
Constraint_BurgersControl::compute_norm
Real compute_norm(const std::vector< Real > &r)
Definition: example_05.hpp:71
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:94
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:311
ROL::Constraint_SimOpt
Defines the constraint operator interface for simulation-based optimization.
Definition: ROL_Constraint_SimOpt.hpp:103
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:192
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:170
ROL_MonteCarloGenerator.hpp
ROL_Objective_SimOpt.hpp
Generated by   doxygen 1.8.5