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