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