ROL
example_08.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Rapid Optimization Library (ROL) Package
4 //
5 // Copyright 2014 NTESS and the ROL contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
15 #include "ROL_Types.hpp"
16 #include "ROL_StdVector.hpp"
17 #include "ROL_BoundConstraint.hpp"
19 #include "ROL_Objective_SimOpt.hpp"
20 #include "ROL_TeuchosBatchManager.hpp"
21 
22 #include "Teuchos_LAPACK.hpp"
23 
24 template<class Real>
25 class L2VectorPrimal;
26 
27 template<class Real>
28 class L2VectorDual;
29 
30 template<class Real>
31 class H1VectorPrimal;
32 
33 template<class Real>
34 class H1VectorDual;
35 
36 template<class Real>
37 class BurgersFEM {
38 private:
39  int nx_;
40  Real dx_;
41  Real nu_;
42  Real nl_;
43  Real u0_;
44  Real u1_;
45  Real f_;
46  Real cH1_;
47  Real cL2_;
48 
49 private:
50  void update(std::vector<Real> &u, const std::vector<Real> &s, const Real alpha=1.0) const {
51  for (unsigned i=0; i<u.size(); i++) {
52  u[i] += alpha*s[i];
53  }
54  }
55 
56  void axpy(std::vector<Real> &out, const Real a, const std::vector<Real> &x, const std::vector<Real> &y) const {
57  for (unsigned i=0; i < x.size(); i++) {
58  out[i] = a*x[i] + y[i];
59  }
60  }
61 
62  void scale(std::vector<Real> &u, const Real alpha=0.0) const {
63  for (unsigned i=0; i<u.size(); i++) {
64  u[i] *= alpha;
65  }
66  }
67 
68  void linear_solve(std::vector<Real> &u, std::vector<Real> &dl, std::vector<Real> &d, std::vector<Real> &du,
69  const std::vector<Real> &r, const bool transpose = false) const {
70  if ( r.size() == 1 ) {
71  u.resize(1,r[0]/d[0]);
72  }
73  else if ( r.size() == 2 ) {
74  u.resize(2,0.0);
75  Real det = d[0]*d[1] - dl[0]*du[0];
76  u[0] = (d[1]*r[0] - du[0]*r[1])/det;
77  u[1] = (d[0]*r[1] - dl[0]*r[0])/det;
78  }
79  else {
80  u.assign(r.begin(),r.end());
81  // Perform LDL factorization
82  Teuchos::LAPACK<int,Real> lp;
83  std::vector<Real> du2(r.size()-2,0.0);
84  std::vector<int> ipiv(r.size(),0);
85  int info;
86  int dim = r.size();
87  int ldb = r.size();
88  int nhrs = 1;
89  lp.GTTRF(dim,&dl[0],&d[0],&du[0],&du2[0],&ipiv[0],&info);
90  char trans = 'N';
91  if ( transpose ) {
92  trans = 'T';
93  }
94  lp.GTTRS(trans,dim,nhrs,&dl[0],&d[0],&du[0],&du2[0],&ipiv[0],&u[0],ldb,&info);
95  }
96  }
97 
98 public:
99  BurgersFEM(int nx = 128, Real nl = 1.0, Real cH1 = 1.0, Real cL2 = 1.0)
100  : nx_(nx), dx_(1.0/((Real)nx+1.0)), nl_(nl), cH1_(cH1), cL2_(cL2) {}
101 
102  void set_problem_data(const Real nu, const Real f, const Real u0, const Real u1) {
103  nu_ = std::pow(10.0,nu-2.0);
104  f_ = 0.01*f;
105  u0_ = 1.0+0.001*u0;
106  u1_ = 0.001*u1;
107  }
108 
109  int num_dof(void) const {
110  return nx_;
111  }
112 
113  Real mesh_spacing(void) const {
114  return dx_;
115  }
116 
117  /***************************************************************************/
118  /*********************** L2 INNER PRODUCT FUNCTIONS ************************/
119  /***************************************************************************/
120  // Compute L2 inner product
121  Real compute_L2_dot(const std::vector<Real> &x, const std::vector<Real> &y) const {
122  Real ip = 0.0;
123  Real c = (((int)x.size()==nx_) ? 4.0 : 2.0);
124  for (unsigned i=0; i<x.size(); i++) {
125  if ( i == 0 ) {
126  ip += dx_/6.0*(c*x[i] + x[i+1])*y[i];
127  }
128  else if ( i == x.size()-1 ) {
129  ip += dx_/6.0*(x[i-1] + c*x[i])*y[i];
130  }
131  else {
132  ip += dx_/6.0*(x[i-1] + 4.0*x[i] + x[i+1])*y[i];
133  }
134  }
135  return ip;
136  }
137 
138  // compute L2 norm
139  Real compute_L2_norm(const std::vector<Real> &r) const {
140  return std::sqrt(compute_L2_dot(r,r));
141  }
142 
143  // Apply L2 Reisz operator
144  void apply_mass(std::vector<Real> &Mu, const std::vector<Real> &u ) const {
145  Mu.resize(u.size(),0.0);
146  Real c = (((int)u.size()==nx_) ? 4.0 : 2.0);
147  for (unsigned i=0; i<u.size(); i++) {
148  if ( i == 0 ) {
149  Mu[i] = dx_/6.0*(c*u[i] + u[i+1]);
150  }
151  else if ( i == u.size()-1 ) {
152  Mu[i] = dx_/6.0*(u[i-1] + c*u[i]);
153  }
154  else {
155  Mu[i] = dx_/6.0*(u[i-1] + 4.0*u[i] + u[i+1]);
156  }
157  }
158  }
159 
160  // Apply L2 inverse Reisz operator
161  void apply_inverse_mass(std::vector<Real> &Mu, const std::vector<Real> &u) const {
162  unsigned nx = u.size();
163  // Build mass matrix
164  std::vector<Real> dl(nx-1,dx_/6.0);
165  std::vector<Real> d(nx,2.0*dx_/3.0);
166  std::vector<Real> du(nx-1,dx_/6.0);
167  if ( (int)nx != nx_ ) {
168  d[ 0] = dx_/3.0;
169  d[nx-1] = dx_/3.0;
170  }
171  linear_solve(Mu,dl,d,du,u);
172  }
173 
174  void test_inverse_mass(std::ostream &outStream = std::cout) {
175  // State Mass Matrix
176  std::vector<Real> u(nx_,0.0), Mu(nx_,0.0), iMMu(nx_,0.0), diff(nx_,0.0);
177  for (int i = 0; i < nx_; i++) {
178  u[i] = 2.0*(Real)rand()/(Real)RAND_MAX - 1.0;
179  }
180  apply_mass(Mu,u);
181  apply_inverse_mass(iMMu,Mu);
182  axpy(diff,-1.0,iMMu,u);
183  Real error = compute_L2_norm(diff);
184  Real normu = compute_L2_norm(u);
185  outStream << "Test Inverse State Mass Matrix\n";
186  outStream << " ||u - inv(M)Mu|| = " << error << "\n";
187  outStream << " ||u|| = " << normu << "\n";
188  outStream << " Relative Error = " << error/normu << "\n";
189  outStream << "\n";
190  // Control Mass Matrix
191  u.resize(nx_+2,0.0); Mu.resize(nx_+2,0.0); iMMu.resize(nx_+2,0.0); diff.resize(nx_+2,0.0);
192  for (int i = 0; i < nx_+2; i++) {
193  u[i] = 2.0*(Real)rand()/(Real)RAND_MAX - 1.0;
194  }
195  apply_mass(Mu,u);
196  apply_inverse_mass(iMMu,Mu);
197  axpy(diff,-1.0,iMMu,u);
198  error = compute_L2_norm(diff);
199  normu = compute_L2_norm(u);
200  outStream << "Test Inverse Control Mass Matrix\n";
201  outStream << " ||z - inv(M)Mz|| = " << error << "\n";
202  outStream << " ||z|| = " << normu << "\n";
203  outStream << " Relative Error = " << error/normu << "\n";
204  outStream << "\n";
205  }
206 
207  /***************************************************************************/
208  /*********************** H1 INNER PRODUCT FUNCTIONS ************************/
209  /***************************************************************************/
210  // Compute H1 inner product
211  Real compute_H1_dot(const std::vector<Real> &x, const std::vector<Real> &y) const {
212  Real ip = 0.0;
213  for (int i=0; i<nx_; i++) {
214  if ( i == 0 ) {
215  ip += cL2_*dx_/6.0*(4.0*x[i] + x[i+1])*y[i]; // Mass term
216  ip += cH1_*(2.0*x[i] - x[i+1])/dx_*y[i]; // Stiffness term
217  }
218  else if ( i == nx_-1 ) {
219  ip += cL2_*dx_/6.0*(x[i-1] + 4.0*x[i])*y[i]; // Mass term
220  ip += cH1_*(2.0*x[i] - x[i-1])/dx_*y[i]; // Stiffness term
221  }
222  else {
223  ip += cL2_*dx_/6.0*(x[i-1] + 4.0*x[i] + x[i+1])*y[i]; // Mass term
224  ip += cH1_*(2.0*x[i] - x[i-1] - x[i+1])/dx_*y[i]; // Stiffness term
225  }
226  }
227  return ip;
228  }
229 
230  // compute H1 norm
231  Real compute_H1_norm(const std::vector<Real> &r) const {
232  return std::sqrt(compute_H1_dot(r,r));
233  }
234 
235  // Apply H2 Reisz operator
236  void apply_H1(std::vector<Real> &Mu, const std::vector<Real> &u ) const {
237  Mu.resize(nx_,0.0);
238  for (int i=0; i<nx_; i++) {
239  if ( i == 0 ) {
240  Mu[i] = cL2_*dx_/6.0*(4.0*u[i] + u[i+1])
241  + cH1_*(2.0*u[i] - u[i+1])/dx_;
242  }
243  else if ( i == nx_-1 ) {
244  Mu[i] = cL2_*dx_/6.0*(u[i-1] + 4.0*u[i])
245  + cH1_*(2.0*u[i] - u[i-1])/dx_;
246  }
247  else {
248  Mu[i] = cL2_*dx_/6.0*(u[i-1] + 4.0*u[i] + u[i+1])
249  + cH1_*(2.0*u[i] - u[i-1] - u[i+1])/dx_;
250  }
251  }
252  }
253 
254  // Apply H1 inverse Reisz operator
255  void apply_inverse_H1(std::vector<Real> &Mu, const std::vector<Real> &u) const {
256  // Build mass matrix
257  std::vector<Real> dl(nx_-1,cL2_*dx_/6.0 - cH1_/dx_);
258  std::vector<Real> d(nx_,2.0*(cL2_*dx_/3.0 + cH1_/dx_));
259  std::vector<Real> du(nx_-1,cL2_*dx_/6.0 - cH1_/dx_);
260  linear_solve(Mu,dl,d,du,u);
261  }
262 
263  void test_inverse_H1(std::ostream &outStream = std::cout) {
264  std::vector<Real> u(nx_,0.0), Mu(nx_,0.0), iMMu(nx_,0.0), diff(nx_,0.0);
265  for (int i = 0; i < nx_; i++) {
266  u[i] = 2.0*(Real)rand()/(Real)RAND_MAX - 1.0;
267  }
268  apply_H1(Mu,u);
269  apply_inverse_H1(iMMu,Mu);
270  axpy(diff,-1.0,iMMu,u);
271  Real error = compute_H1_norm(diff);
272  Real normu = compute_H1_norm(u);
273  outStream << "Test Inverse State H1 Matrix\n";
274  outStream << " ||u - inv(M)Mu|| = " << error << "\n";
275  outStream << " ||u|| = " << normu << "\n";
276  outStream << " Relative Error = " << error/normu << "\n";
277  outStream << "\n";
278  }
279 
280  /***************************************************************************/
281  /*********************** PDE RESIDUAL AND SOLVE ****************************/
282  /***************************************************************************/
283  // Compute current PDE residual
284  void compute_residual(std::vector<Real> &r, const std::vector<Real> &u,
285  const std::vector<Real> &z) const {
286  r.clear();
287  r.resize(nx_,0.0);
288  for (int i=0; i<nx_; i++) {
289  // Contribution from stiffness term
290  if ( i==0 ) {
291  r[i] = nu_/dx_*(2.0*u[i]-u[i+1]);
292  }
293  else if (i==nx_-1) {
294  r[i] = nu_/dx_*(2.0*u[i]-u[i-1]);
295  }
296  else {
297  r[i] = nu_/dx_*(2.0*u[i]-u[i-1]-u[i+1]);
298  }
299  // Contribution from nonlinear term
300  if (i<nx_-1){
301  r[i] += nl_*u[i+1]*(u[i]+u[i+1])/6.0;
302  }
303  if (i>0) {
304  r[i] -= nl_*u[i-1]*(u[i-1]+u[i])/6.0;
305  }
306  // Contribution from control
307  r[i] -= dx_/6.0*(z[i]+4.0*z[i+1]+z[i+2]);
308  // Contribution from right-hand side
309  r[i] -= dx_*f_;
310  }
311  // Contribution from Dirichlet boundary terms
312  r[0] -= nl_*(u0_*u[ 0]/6.0 + u0_*u0_/6.0) + nu_*u0_/dx_;
313  r[nx_-1] += nl_*(u1_*u[nx_-1]/6.0 + u1_*u1_/6.0) - nu_*u1_/dx_;
314  }
315 
316  /***************************************************************************/
317  /*********************** PDE JACOBIAN FUNCTIONS ****************************/
318  /***************************************************************************/
319  // Build PDE Jacobian trigiagonal matrix
320  void compute_pde_jacobian(std::vector<Real> &dl, // Lower diagonal
321  std::vector<Real> &d, // Diagonal
322  std::vector<Real> &du, // Upper diagonal
323  const std::vector<Real> &u) const { // State variable
324  // Get Diagonal and Off-Diagonal Entries of linear PDE Jacobian
325  d.clear();
326  d.resize(nx_,nu_*2.0/dx_);
327  dl.clear();
328  dl.resize(nx_-1,-nu_/dx_);
329  du.clear();
330  du.resize(nx_-1,-nu_/dx_);
331  // Contribution from nonlinearity
332  for (int i=0; i<nx_; i++) {
333  if (i<nx_-1) {
334  dl[i] += nl_*(-2.0*u[i]-u[i+1])/6.0;
335  d[i] += nl_*u[i+1]/6.0;
336  }
337  if (i>0) {
338  d[i] -= nl_*u[i-1]/6.0;
339  du[i-1] += nl_*(u[i-1]+2.0*u[i])/6.0;
340  }
341  }
342  // Contribution from Dirichlet boundary conditions
343  d[ 0] -= nl_*u0_/6.0;
344  d[nx_-1] += nl_*u1_/6.0;
345  }
346 
347  // Apply PDE Jacobian to a vector
348  void apply_pde_jacobian(std::vector<Real> &jv,
349  const std::vector<Real> &v,
350  const std::vector<Real> &u,
351  const std::vector<Real> &z) const {
352  // Fill jv
353  for (int i = 0; i < nx_; i++) {
354  jv[i] = nu_/dx_*2.0*v[i];
355  if ( i > 0 ) {
356  jv[i] += -nu_/dx_*v[i-1]-nl_*(u[i-1]/6.0*v[i]+(u[i]+2.0*u[i-1])/6.0*v[i-1]);
357  }
358  if ( i < nx_-1 ) {
359  jv[i] += -nu_/dx_*v[i+1]+nl_*(u[i+1]/6.0*v[i]+(u[i]+2.0*u[i+1])/6.0*v[i+1]);
360  }
361  }
362  jv[ 0] -= nl_*u0_/6.0*v[0];
363  jv[nx_-1] += nl_*u1_/6.0*v[nx_-1];
364  }
365 
366  // Apply inverse PDE Jacobian to a vector
367  void apply_inverse_pde_jacobian(std::vector<Real> &ijv,
368  const std::vector<Real> &v,
369  const std::vector<Real> &u,
370  const std::vector<Real> &z) const {
371  // Get PDE Jacobian
372  std::vector<Real> d(nx_,0.0);
373  std::vector<Real> dl(nx_-1,0.0);
374  std::vector<Real> du(nx_-1,0.0);
375  compute_pde_jacobian(dl,d,du,u);
376  // Solve solve state sensitivity system at current time step
377  linear_solve(ijv,dl,d,du,v);
378  }
379 
380  // Apply adjoint PDE Jacobian to a vector
381  void apply_adjoint_pde_jacobian(std::vector<Real> &ajv,
382  const std::vector<Real> &v,
383  const std::vector<Real> &u,
384  const std::vector<Real> &z) const {
385  // Fill jvp
386  for (int i = 0; i < nx_; i++) {
387  ajv[i] = nu_/dx_*2.0*v[i];
388  if ( i > 0 ) {
389  ajv[i] += -nu_/dx_*v[i-1]-nl_*(u[i-1]/6.0*v[i]
390  -(u[i-1]+2.0*u[i])/6.0*v[i-1]);
391  }
392  if ( i < nx_-1 ) {
393  ajv[i] += -nu_/dx_*v[i+1]+nl_*(u[i+1]/6.0*v[i]
394  -(u[i+1]+2.0*u[i])/6.0*v[i+1]);
395  }
396  }
397  ajv[ 0] -= nl_*u0_/6.0*v[0];
398  ajv[nx_-1] += nl_*u1_/6.0*v[nx_-1];
399  }
400 
401  // Apply inverse adjoint PDE Jacobian to a vector
402  void apply_inverse_adjoint_pde_jacobian(std::vector<Real> &iajv,
403  const std::vector<Real> &v,
404  const std::vector<Real> &u,
405  const std::vector<Real> &z) const {
406  // Get PDE Jacobian
407  std::vector<Real> d(nx_,0.0);
408  std::vector<Real> du(nx_-1,0.0);
409  std::vector<Real> dl(nx_-1,0.0);
410  compute_pde_jacobian(dl,d,du,u);
411  // Solve solve adjoint system at current time step
412  linear_solve(iajv,dl,d,du,v,true);
413  }
414 
415  /***************************************************************************/
416  /*********************** CONTROL JACOBIAN FUNCTIONS ************************/
417  /***************************************************************************/
418  // Apply control Jacobian to a vector
419  void apply_control_jacobian(std::vector<Real> &jv,
420  const std::vector<Real> &v,
421  const std::vector<Real> &u,
422  const std::vector<Real> &z) const {
423  for (int i=0; i<nx_; i++) {
424  // Contribution from control
425  jv[i] = -dx_/6.0*(v[i]+4.0*v[i+1]+v[i+2]);
426  }
427  }
428 
429  // Apply adjoint control Jacobian to a vector
430  void apply_adjoint_control_jacobian(std::vector<Real> &jv,
431  const std::vector<Real> &v,
432  const std::vector<Real> &u,
433  const std::vector<Real> &z) const {
434  for (int i=0; i<nx_+2; i++) {
435  if ( i == 0 ) {
436  jv[i] = -dx_/6.0*v[i];
437  }
438  else if ( i == 1 ) {
439  jv[i] = -dx_/6.0*(4.0*v[i-1]+v[i]);
440  }
441  else if ( i == nx_ ) {
442  jv[i] = -dx_/6.0*(4.0*v[i-1]+v[i-2]);
443  }
444  else if ( i == nx_+1 ) {
445  jv[i] = -dx_/6.0*v[i-2];
446  }
447  else {
448  jv[i] = -dx_/6.0*(v[i-2]+4.0*v[i-1]+v[i]);
449  }
450  }
451  }
452 
453  /***************************************************************************/
454  /*********************** AJDOINT HESSIANS **********************************/
455  /***************************************************************************/
456  void apply_adjoint_pde_hessian(std::vector<Real> &ahwv,
457  const std::vector<Real> &w,
458  const std::vector<Real> &v,
459  const std::vector<Real> &u,
460  const std::vector<Real> &z) const {
461  for (int i=0; i<nx_; i++) {
462  // Contribution from nonlinear term
463  ahwv[i] = 0.0;
464  if (i<nx_-1){
465  ahwv[i] += (w[i]*v[i+1] - w[i+1]*(2.0*v[i]+v[i+1]))/6.0;
466  }
467  if (i>0) {
468  ahwv[i] += (w[i-1]*(v[i-1]+2.0*v[i]) - w[i]*v[i-1])/6.0;
469  }
470  }
471  //ahwv.assign(u.size(),0.0);
472  }
473  void apply_adjoint_pde_control_hessian(std::vector<Real> &ahwv,
474  const std::vector<Real> &w,
475  const std::vector<Real> &v,
476  const std::vector<Real> &u,
477  const std::vector<Real> &z) {
478  ahwv.assign(u.size(),0.0);
479  }
480  void apply_adjoint_control_pde_hessian(std::vector<Real> &ahwv,
481  const std::vector<Real> &w,
482  const std::vector<Real> &v,
483  const std::vector<Real> &u,
484  const std::vector<Real> &z) {
485  ahwv.assign(z.size(),0.0);
486  }
487  void apply_adjoint_control_hessian(std::vector<Real> &ahwv,
488  const std::vector<Real> &w,
489  const std::vector<Real> &v,
490  const std::vector<Real> &u,
491  const std::vector<Real> &z) {
492  ahwv.assign(z.size(),0.0);
493  }
494 };
495 
496 template<class Real>
497 class L2VectorPrimal : public ROL::StdVector<Real> {
498 private:
499  ROL::Ptr<BurgersFEM<Real> > fem_;
500  mutable ROL::Ptr<L2VectorDual<Real> > dual_vec_;
501  mutable bool isDualInitialized_;
502 
503 public:
504  L2VectorPrimal(const ROL::Ptr<std::vector<Real> > & vec,
505  const ROL::Ptr<BurgersFEM<Real> > &fem)
506  : ROL::StdVector<Real>(vec),
507  fem_(fem), dual_vec_(ROL::nullPtr), isDualInitialized_(false) {}
508 
509  Real dot( const ROL::Vector<Real> &x ) const {
510  const L2VectorPrimal & ex = dynamic_cast<const L2VectorPrimal&>(x);
511  const std::vector<Real>& xval = *ex.getVector();
512  const std::vector<Real>& yval = *(ROL::StdVector<Real>::getVector());
513  return fem_->compute_L2_dot(xval,yval);
514  }
515 
516  ROL::Ptr<ROL::Vector<Real> > clone() const {
517  return ROL::makePtr<L2VectorPrimal>(
518  ROL::makePtr<std::vector<Real>>(ROL::StdVector<Real>::dimension()),fem_);
519  }
520 
521  const ROL::Vector<Real>& dual() const {
522  if ( !isDualInitialized_) {
523  dual_vec_ = ROL::makePtr<L2VectorDual<Real>>(
524  ROL::makePtr<std::vector<Real>>(ROL::StdVector<Real>::dimension()),fem_);
525  isDualInitialized_ = true;
526  }
527  fem_->apply_mass(*(dual_vec_->getVector()),*(ROL::StdVector<Real>::getVector()));
528  return *dual_vec_;
529  }
530 
531 };
532 
533 template<class Real>
534 class L2VectorDual : public ROL::StdVector<Real> {
535 private:
536  ROL::Ptr<BurgersFEM<Real> > fem_;
537  mutable ROL::Ptr<L2VectorPrimal<Real> > prim_vec_;
538  mutable bool isPrimalInitialized_;
539 
540 public:
541  L2VectorDual(const ROL::Ptr<std::vector<Real> > & vec,
542  const ROL::Ptr<BurgersFEM<Real> > &fem)
543  : ROL::StdVector<Real>(vec),
544  fem_(fem), prim_vec_(ROL::nullPtr), isPrimalInitialized_(false) {}
545 
546  Real dot( const ROL::Vector<Real> &x ) const {
547  const L2VectorDual & ex = dynamic_cast<const L2VectorDual&>(x);
548  const std::vector<Real>& xval = *ex.getVector();
549  const std::vector<Real>& yval = *(ROL::StdVector<Real>::getVector());
550  unsigned dimension = yval.size();
551  std::vector<Real> Mx(dimension,0.0);
552  fem_->apply_inverse_mass(Mx,xval);
553  Real val(0);
554  for (unsigned i = 0; i < dimension; i++) {
555  val += Mx[i]*yval[i];
556  }
557  return val;
558  }
559 
560  ROL::Ptr<ROL::Vector<Real> > clone() const {
561  return ROL::makePtr<L2VectorDual>(
562  ROL::makePtr<std::vector<Real>>(ROL::StdVector<Real>::dimension()),fem_);
563  }
564 
565  const ROL::Vector<Real>& dual() const {
566  if (!isPrimalInitialized_) {
567  prim_vec_ = ROL::makePtr<L2VectorPrimal<Real>>(
568  ROL::makePtr<std::vector<Real>>(ROL::StdVector<Real>::dimension()),fem_);
569  isPrimalInitialized_ = true;
570  }
571  fem_->apply_inverse_mass(*(prim_vec_->getVector()),*(ROL::StdVector<Real>::getVector()));
572  return *prim_vec_;
573  }
574 
575 };
576 
577 template<class Real>
578 class H1VectorPrimal : public ROL::StdVector<Real> {
579 private:
580  ROL::Ptr<BurgersFEM<Real> > fem_;
581  mutable ROL::Ptr<H1VectorDual<Real> > dual_vec_;
582  mutable bool isDualInitialized_;
583 
584 public:
585  H1VectorPrimal(const ROL::Ptr<std::vector<Real> > & vec,
586  const ROL::Ptr<BurgersFEM<Real> > &fem)
587  : ROL::StdVector<Real>(vec),
588  fem_(fem), dual_vec_(ROL::nullPtr), isDualInitialized_(false) {}
589 
590  Real dot( const ROL::Vector<Real> &x ) const {
591  const H1VectorPrimal & ex = dynamic_cast<const H1VectorPrimal&>(x);
592  const std::vector<Real>& xval = *ex.getVector();
593  const std::vector<Real>& yval = *(ROL::StdVector<Real>::getVector());
594  return fem_->compute_H1_dot(xval,yval);
595  }
596 
597  ROL::Ptr<ROL::Vector<Real> > clone() const {
598  return ROL::makePtr<H1VectorPrimal>(
599  ROL::makePtr<std::vector<Real>>(ROL::StdVector<Real>::dimension()),fem_);
600  }
601 
602  const ROL::Vector<Real>& dual() const {
603  if ( !isDualInitialized_) {
604  dual_vec_ = ROL::makePtr<H1VectorDual<Real>>(
605  ROL::makePtr<std::vector<Real>>(ROL::StdVector<Real>::dimension()),fem_);
606  isDualInitialized_ = true;
607  }
608  fem_->apply_H1(*(dual_vec_->getVector()),*(ROL::StdVector<Real>::getVector()));
609  return *dual_vec_;
610  }
611 
612 };
613 
614 template<class Real>
615 class H1VectorDual : public ROL::StdVector<Real> {
616 private:
617  ROL::Ptr<BurgersFEM<Real> > fem_;
618  mutable ROL::Ptr<H1VectorPrimal<Real> > prim_vec_;
619  mutable bool isPrimalInitialized_;
620 
621 public:
622  H1VectorDual(const ROL::Ptr<std::vector<Real> > & vec,
623  const ROL::Ptr<BurgersFEM<Real> > &fem)
624  : ROL::StdVector<Real>(vec),
625  fem_(fem), prim_vec_(ROL::nullPtr), isPrimalInitialized_(false) {}
626 
627  Real dot( const ROL::Vector<Real> &x ) const {
628  const H1VectorDual & ex = dynamic_cast<const H1VectorDual&>(x);
629  const std::vector<Real>& xval = *ex.getVector();
630  const std::vector<Real>& yval = *(ROL::StdVector<Real>::getVector());
631  unsigned dimension = yval.size();
632  std::vector<Real> Mx(dimension,0.0);
633  fem_->apply_inverse_H1(Mx,xval);
634  Real val(0);
635  for (unsigned i = 0; i < dimension; i++) {
636  val += Mx[i]*yval[i];
637  }
638  return val;
639  }
640 
641  ROL::Ptr<ROL::Vector<Real> > clone() const {
642  return ROL::makePtr<H1VectorDual>(
643  ROL::makePtr<std::vector<Real>>(ROL::StdVector<Real>::dimension()),fem_);
644  }
645 
646  const ROL::Vector<Real>& dual() const {
647  if (!isPrimalInitialized_) {
648  prim_vec_ = ROL::makePtr<H1VectorPrimal<Real>>(
649  ROL::makePtr<std::vector<Real>>(ROL::StdVector<Real>::dimension()),fem_);
650  isPrimalInitialized_ = true;
651  }
652  fem_->apply_inverse_H1(*(prim_vec_->getVector()),*(ROL::StdVector<Real>::getVector()));
653  return *prim_vec_;
654  }
655 
656 };
657 
658 template<class Real>
660 private:
661 
664 
667 
670 
671  ROL::Ptr<BurgersFEM<Real> > fem_;
672  bool useHessian_;
673 
674 public:
675  Constraint_BurgersControl(ROL::Ptr<BurgersFEM<Real> > &fem, bool useHessian = true)
676  : fem_(fem), useHessian_(useHessian) {}
677 
679  const ROL::Vector<Real> &z, Real &tol) {
680  ROL::Ptr<std::vector<Real> > cp =
681  dynamic_cast<PrimalConstraintVector&>(c).getVector();
682  ROL::Ptr<const std::vector<Real> > up =
683  dynamic_cast<const PrimalStateVector&>(u).getVector();
684  ROL::Ptr<const std::vector<Real> > zp =
685  dynamic_cast<const PrimalControlVector&>(z).getVector();
686 
687  const std::vector<Real> param
689  fem_->set_problem_data(param[0],param[1],param[2],param[3]);
690 
691  fem_->compute_residual(*cp,*up,*zp);
692  }
693 
695  const ROL::Vector<Real> &z, Real &tol) {
696  ROL::Ptr<std::vector<Real> > jvp =
697  dynamic_cast<PrimalConstraintVector&>(jv).getVector();
698  ROL::Ptr<const std::vector<Real> > vp =
699  dynamic_cast<const PrimalStateVector&>(v).getVector();
700  ROL::Ptr<const std::vector<Real> > up =
701  dynamic_cast<const PrimalStateVector&>(u).getVector();
702  ROL::Ptr<const std::vector<Real> > zp =
703  dynamic_cast<const PrimalControlVector&>(z).getVector();
704 
705  const std::vector<Real> param
707  fem_->set_problem_data(param[0],param[1],param[2],param[3]);
708 
709  fem_->apply_pde_jacobian(*jvp,*vp,*up,*zp);
710  }
711 
713  const ROL::Vector<Real> &z, Real &tol) {
714  ROL::Ptr<std::vector<Real> > jvp =
715  dynamic_cast<PrimalConstraintVector&>(jv).getVector();
716  ROL::Ptr<const std::vector<Real> > vp =
717  dynamic_cast<const PrimalControlVector&>(v).getVector();
718  ROL::Ptr<const std::vector<Real> > up =
719  dynamic_cast<const PrimalStateVector&>(u).getVector();
720  ROL::Ptr<const std::vector<Real> > zp =
721  dynamic_cast<const PrimalControlVector&>(z).getVector();
722 
723  const std::vector<Real> param
725  fem_->set_problem_data(param[0],param[1],param[2],param[3]);
726 
727  fem_->apply_control_jacobian(*jvp,*vp,*up,*zp);
728  }
729 
731  const ROL::Vector<Real> &z, Real &tol) {
732  ROL::Ptr<std::vector<Real> > ijvp =
733  dynamic_cast<PrimalStateVector&>(ijv).getVector();
734  ROL::Ptr<const std::vector<Real> > vp =
735  dynamic_cast<const PrimalConstraintVector&>(v).getVector();
736  ROL::Ptr<const std::vector<Real> > up =
737  dynamic_cast<const PrimalStateVector&>(u).getVector();
738  ROL::Ptr<const std::vector<Real> > zp =
739  dynamic_cast<const PrimalControlVector&>(z).getVector();
740 
741  const std::vector<Real> param
743  fem_->set_problem_data(param[0],param[1],param[2],param[3]);
744 
745  fem_->apply_inverse_pde_jacobian(*ijvp,*vp,*up,*zp);
746  }
747 
749  const ROL::Vector<Real> &z, Real &tol) {
750  ROL::Ptr<std::vector<Real> > jvp =
751  dynamic_cast<DualStateVector&>(ajv).getVector();
752  ROL::Ptr<const std::vector<Real> > vp =
753  dynamic_cast<const DualConstraintVector&>(v).getVector();
754  ROL::Ptr<const std::vector<Real> > up =
755  dynamic_cast<const PrimalStateVector&>(u).getVector();
756  ROL::Ptr<const std::vector<Real> > zp =
757  dynamic_cast<const PrimalControlVector&>(z).getVector();
758 
759  const std::vector<Real> param
761  fem_->set_problem_data(param[0],param[1],param[2],param[3]);
762 
763  fem_->apply_adjoint_pde_jacobian(*jvp,*vp,*up,*zp);
764  }
765 
767  const ROL::Vector<Real> &z, Real &tol) {
768  ROL::Ptr<std::vector<Real> > jvp =
769  dynamic_cast<DualControlVector&>(jv).getVector();
770  ROL::Ptr<const std::vector<Real> > vp =
771  dynamic_cast<const DualConstraintVector&>(v).getVector();
772  ROL::Ptr<const std::vector<Real> > up =
773  dynamic_cast<const PrimalStateVector&>(u).getVector();
774  ROL::Ptr<const std::vector<Real> > zp =
775  dynamic_cast<const PrimalControlVector&>(z).getVector();
776 
777  const std::vector<Real> param
779  fem_->set_problem_data(param[0],param[1],param[2],param[3]);
780 
781  fem_->apply_adjoint_control_jacobian(*jvp,*vp,*up,*zp);
782  }
783 
785  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
786  ROL::Ptr<std::vector<Real> > iajvp =
787  dynamic_cast<DualConstraintVector&>(iajv).getVector();
788  ROL::Ptr<const std::vector<Real> > vp =
789  dynamic_cast<const DualStateVector&>(v).getVector();
790  ROL::Ptr<const std::vector<Real> > up =
791  dynamic_cast<const PrimalStateVector&>(u).getVector();
792  ROL::Ptr<const std::vector<Real> > zp =
793  dynamic_cast<const PrimalControlVector&>(z).getVector();
794 
795  const std::vector<Real> param
797  fem_->set_problem_data(param[0],param[1],param[2],param[3]);
798 
799  fem_->apply_inverse_adjoint_pde_jacobian(*iajvp,*vp,*up,*zp);
800  }
801 
803  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
804  if ( useHessian_ ) {
805  ROL::Ptr<std::vector<Real> > ahwvp =
806  dynamic_cast<DualStateVector&>(ahwv).getVector();
807  ROL::Ptr<const std::vector<Real> > wp =
808  dynamic_cast<const DualConstraintVector&>(w).getVector();
809  ROL::Ptr<const std::vector<Real> > vp =
810  dynamic_cast<const PrimalStateVector&>(v).getVector();
811  ROL::Ptr<const std::vector<Real> > up =
812  dynamic_cast<const PrimalStateVector&>(u).getVector();
813  ROL::Ptr<const std::vector<Real> > zp =
814  dynamic_cast<const PrimalControlVector&>(z).getVector();
815 
816  const std::vector<Real> param
818  fem_->set_problem_data(param[0],param[1],param[2],param[3]);
819 
820  fem_->apply_adjoint_pde_hessian(*ahwvp,*wp,*vp,*up,*zp);
821  }
822  else {
823  ahwv.zero();
824  }
825  }
826 
828  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
829  if ( useHessian_ ) {
830  ROL::Ptr<std::vector<Real> > ahwvp =
831  dynamic_cast<DualControlVector&>(ahwv).getVector();
832  ROL::Ptr<const std::vector<Real> > wp =
833  dynamic_cast<const DualConstraintVector&>(w).getVector();
834  ROL::Ptr<const std::vector<Real> > vp =
835  dynamic_cast<const PrimalStateVector&>(v).getVector();
836  ROL::Ptr<const std::vector<Real> > up =
837  dynamic_cast<const PrimalStateVector&>(u).getVector();
838  ROL::Ptr<const std::vector<Real> > zp =
839  dynamic_cast<const PrimalControlVector&>(z).getVector();
840 
841  const std::vector<Real> param
843  fem_->set_problem_data(param[0],param[1],param[2],param[3]);
844 
845  fem_->apply_adjoint_control_pde_hessian(*ahwvp,*wp,*vp,*up,*zp);
846  }
847  else {
848  ahwv.zero();
849  }
850  }
852  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
853  if ( useHessian_ ) {
854  ROL::Ptr<std::vector<Real> > ahwvp =
855  dynamic_cast<DualStateVector&>(ahwv).getVector();
856  ROL::Ptr<const std::vector<Real> > wp =
857  dynamic_cast<const DualConstraintVector&>(w).getVector();
858  ROL::Ptr<const std::vector<Real> > vp =
859  dynamic_cast<const PrimalControlVector&>(v).getVector();
860  ROL::Ptr<const std::vector<Real> > up =
861  dynamic_cast<const PrimalStateVector&>(u).getVector();
862  ROL::Ptr<const std::vector<Real> > zp =
863  dynamic_cast<const PrimalControlVector&>(z).getVector();
864 
865  const std::vector<Real> param
867  fem_->set_problem_data(param[0],param[1],param[2],param[3]);
868 
869  fem_->apply_adjoint_pde_control_hessian(*ahwvp,*wp,*vp,*up,*zp);
870  }
871  else {
872  ahwv.zero();
873  }
874  }
876  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
877  if ( useHessian_ ) {
878  ROL::Ptr<std::vector<Real> > ahwvp =
879  dynamic_cast<DualControlVector&>(ahwv).getVector();
880  ROL::Ptr<const std::vector<Real> > wp =
881  dynamic_cast<const DualConstraintVector&>(w).getVector();
882  ROL::Ptr<const std::vector<Real> > vp =
883  dynamic_cast<const PrimalControlVector&>(v).getVector();
884  ROL::Ptr<const std::vector<Real> > up =
885  dynamic_cast<const PrimalStateVector&>(u).getVector();
886  ROL::Ptr<const std::vector<Real> > zp =
887  dynamic_cast<const PrimalControlVector&>(z).getVector();
888 
889  const std::vector<Real> param
891  fem_->set_problem_data(param[0],param[1],param[2],param[3]);
892 
893  fem_->apply_adjoint_control_hessian(*ahwvp,*wp,*vp,*up,*zp);
894  }
895  else {
896  ahwv.zero();
897  }
898  }
899 };
900 
901 template<class Real>
902 class Objective_BurgersControl : public ROL::Objective_SimOpt<Real> {
903 private:
904 
907 
910 
911  Real alpha_; // Penalty Parameter
912  ROL::Ptr<BurgersFEM<Real> > fem_;
913  ROL::Ptr<ROL::Vector<Real> > ud_;
914  ROL::Ptr<ROL::Vector<Real> > diff_;
915 
916 public:
918  const ROL::Ptr<ROL::Vector<Real> > &ud,
919  Real alpha = 1.e-4) : alpha_(alpha), fem_(fem), ud_(ud) {
920  diff_ = ud_->clone();
921  }
922 
924 
925  Real value( const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
926  ROL::Ptr<const std::vector<Real> > up =
927  dynamic_cast<const PrimalStateVector&>(u).getVector();
928  ROL::Ptr<const std::vector<Real> > zp =
929  dynamic_cast<const PrimalControlVector&>(z).getVector();
930  ROL::Ptr<const std::vector<Real> > udp =
931  dynamic_cast<const L2VectorPrimal<Real>& >(*ud_).getVector();
932 
933  std::vector<Real> diff(udp->size(),0.0);
934  for (unsigned i = 0; i < udp->size(); i++) {
935  diff[i] = (*up)[i] - (*udp)[i];
936  }
937  return 0.5*(fem_->compute_L2_dot(diff,diff) + alpha_*fem_->compute_L2_dot(*zp,*zp));
938  }
939 
940  void gradient_1( ROL::Vector<Real> &g, const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
941  ROL::Ptr<std::vector<Real> > gp =
942  dynamic_cast<DualStateVector&>(g).getVector();
943  ROL::Ptr<const std::vector<Real> > up =
944  dynamic_cast<const PrimalStateVector&>(u).getVector();
945  ROL::Ptr<const std::vector<Real> > udp =
946  dynamic_cast<const L2VectorPrimal<Real>&>(*ud_).getVector();
947 
948  std::vector<Real> diff(udp->size(),0.0);
949  for (unsigned i = 0; i < udp->size(); i++) {
950  diff[i] = (*up)[i] - (*udp)[i];
951  }
952  fem_->apply_mass(*gp,diff);
953  }
954 
955  void gradient_2( ROL::Vector<Real> &g, const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
956  ROL::Ptr<std::vector<Real> > gp =
957  dynamic_cast<DualControlVector&>(g).getVector();
958  ROL::Ptr<const std::vector<Real> > zp =
959  dynamic_cast<const PrimalControlVector&>(z).getVector();
960 
961  fem_->apply_mass(*gp,*zp);
962  for (unsigned i = 0; i < zp->size(); i++) {
963  (*gp)[i] *= alpha_;
964  }
965  }
966 
968  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
969  ROL::Ptr<std::vector<Real> > hvp =
970  dynamic_cast<DualStateVector&>(hv).getVector();
971  ROL::Ptr<const std::vector<Real> > vp =
972  dynamic_cast<const PrimalStateVector&>(v).getVector();
973 
974  fem_->apply_mass(*hvp,*vp);
975  }
976 
978  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
979  hv.zero();
980  }
981 
983  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
984  hv.zero();
985  }
986 
988  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
989  ROL::Ptr<std::vector<Real> > hvp =
990  dynamic_cast<DualControlVector&>(hv).getVector();
991  ROL::Ptr<const std::vector<Real> > vp =
992  dynamic_cast<const PrimalControlVector&>(v).getVector();
993 
994  fem_->apply_mass(*hvp,*vp);
995  for (unsigned i = 0; i < vp->size(); i++) {
996  (*hvp)[i] *= alpha_;
997  }
998  }
999 };
1000 
1001 template<class Real>
1002 class L2BoundConstraint : public ROL::BoundConstraint<Real> {
1003 private:
1004  int dim_;
1005  std::vector<Real> x_lo_;
1006  std::vector<Real> x_up_;
1007  Real min_diff_;
1008  Real scale_;
1009  ROL::Ptr<BurgersFEM<Real> > fem_;
1010  ROL::Ptr<ROL::Vector<Real> > l_;
1011  ROL::Ptr<ROL::Vector<Real> > u_;
1012 
1013  void cast_vector(ROL::Ptr<std::vector<Real> > &xvec,
1014  ROL::Vector<Real> &x) const {
1015  try {
1016  xvec = dynamic_cast<L2VectorPrimal<Real>&>(x).getVector();
1017  }
1018  catch (std::exception &e) {
1019  xvec = dynamic_cast<L2VectorDual<Real>&>(x).getVector();
1020  }
1021  }
1022 
1023  void cast_const_vector(ROL::Ptr<const std::vector<Real> > &xvec,
1024  const ROL::Vector<Real> &x) const {
1025  try {
1026  xvec = dynamic_cast<const L2VectorPrimal<Real>&>(x).getVector();
1027  }
1028  catch (std::exception &e) {
1029  xvec = dynamic_cast<const L2VectorDual<Real>&>(x).getVector();
1030  }
1031  }
1032 
1033  void axpy(std::vector<Real> &out, const Real a,
1034  const std::vector<Real> &x, const std::vector<Real> &y) const{
1035  out.resize(dim_,0.0);
1036  for (unsigned i = 0; i < dim_; i++) {
1037  out[i] = a*x[i] + y[i];
1038  }
1039  }
1040 
1041  void projection(std::vector<Real> &x) {
1042  for ( int i = 0; i < dim_; i++ ) {
1043  x[i] = std::max(x_lo_[i],std::min(x_up_[i],x[i]));
1044  }
1045  }
1046 
1047 public:
1048  L2BoundConstraint(std::vector<Real> &l, std::vector<Real> &u,
1049  const ROL::Ptr<BurgersFEM<Real> > &fem, Real scale = 1.0)
1050  : x_lo_(l), x_up_(u), scale_(scale), fem_(fem) {
1051  dim_ = x_lo_.size();
1052  for ( int i = 0; i < dim_; i++ ) {
1053  if ( i == 0 ) {
1054  min_diff_ = x_up_[i] - x_lo_[i];
1055  }
1056  else {
1057  min_diff_ = ( (min_diff_ < (x_up_[i] - x_lo_[i])) ? min_diff_ : (x_up_[i] - x_lo_[i]) );
1058  }
1059  }
1060  min_diff_ *= 0.5;
1061  l_ = ROL::makePtr<L2VectorPrimal<Real>>(
1062  ROL::makePtr<std::vector<Real>>(l), fem);
1063  u_ = ROL::makePtr<L2VectorPrimal<Real>>(
1064  ROL::makePtr<std::vector<Real>>(u), fem);
1065  }
1066 
1067  bool isFeasible( const ROL::Vector<Real> &x ) {
1068  ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1069  bool val = true;
1070  int cnt = 1;
1071  for ( int i = 0; i < dim_; i++ ) {
1072  if ( (*ex)[i] >= x_lo_[i] && (*ex)[i] <= x_up_[i] ) { cnt *= 1; }
1073  else { cnt *= 0; }
1074  }
1075  if ( cnt == 0 ) { val = false; }
1076  return val;
1077  }
1078 
1080  ROL::Ptr<std::vector<Real> > ex; cast_vector(ex,x);
1081  projection(*ex);
1082  }
1083 
1084  void pruneLowerActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &x, Real eps = Real(0)) {
1085  ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1086  ROL::Ptr<std::vector<Real> > ev; cast_vector(ev,v);
1087  Real epsn = std::min(scale_*eps,min_diff_);
1088  for ( int i = 0; i < dim_; i++ ) {
1089  if ( ((*ex)[i] <= x_lo_[i]+epsn) ) {
1090  (*ev)[i] = 0.0;
1091  }
1092  }
1093  }
1094 
1095  void pruneUpperActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &x, Real eps = Real(0)) {
1096  ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1097  ROL::Ptr<std::vector<Real> > ev; cast_vector(ev,v);
1098  Real epsn = std::min(scale_*eps,min_diff_);
1099  for ( int i = 0; i < dim_; i++ ) {
1100  if ( ((*ex)[i] >= x_up_[i]-epsn) ) {
1101  (*ev)[i] = 0.0;
1102  }
1103  }
1104  }
1105 
1106  void pruneActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &x, Real eps = Real(0)) {
1107  ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1108  ROL::Ptr<std::vector<Real> > ev; cast_vector(ev,v);
1109  Real epsn = std::min(scale_*eps,min_diff_);
1110  for ( int i = 0; i < dim_; i++ ) {
1111  if ( ((*ex)[i] <= x_lo_[i]+epsn) ||
1112  ((*ex)[i] >= x_up_[i]-epsn) ) {
1113  (*ev)[i] = 0.0;
1114  }
1115  }
1116  }
1117 
1118  void pruneLowerActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &g, const ROL::Vector<Real> &x, Real xeps = Real(0), Real geps = Real(0)) {
1119  ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1120  ROL::Ptr<const std::vector<Real> > eg; cast_const_vector(eg,g);
1121  ROL::Ptr<std::vector<Real> > ev; cast_vector(ev,v);
1122  Real epsn = std::min(scale_*xeps,min_diff_);
1123  for ( int i = 0; i < dim_; i++ ) {
1124  if ( ((*ex)[i] <= x_lo_[i]+epsn && (*eg)[i] > geps) ) {
1125  (*ev)[i] = 0.0;
1126  }
1127  }
1128  }
1129 
1130  void pruneUpperActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &g, const ROL::Vector<Real> &x, Real xeps = Real(0), Real geps = Real(0)) {
1131  ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1132  ROL::Ptr<const std::vector<Real> > eg; cast_const_vector(eg,g);
1133  ROL::Ptr<std::vector<Real> > ev; cast_vector(ev,v);
1134  Real epsn = std::min(scale_*xeps,min_diff_);
1135  for ( int i = 0; i < dim_; i++ ) {
1136  if ( ((*ex)[i] >= x_up_[i]-epsn && (*eg)[i] < -geps) ) {
1137  (*ev)[i] = 0.0;
1138  }
1139  }
1140  }
1141 
1142  void pruneActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &g, const ROL::Vector<Real> &x, Real xeps = Real(0), Real geps = Real(0)) {
1143  ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1144  ROL::Ptr<const std::vector<Real> > eg; cast_const_vector(eg,g);
1145  ROL::Ptr<std::vector<Real> > ev; cast_vector(ev,v);
1146  Real epsn = std::min(scale_*xeps,min_diff_);
1147  for ( int i = 0; i < dim_; i++ ) {
1148  if ( ((*ex)[i] <= x_lo_[i]+epsn && (*eg)[i] > geps) ||
1149  ((*ex)[i] >= x_up_[i]-epsn && (*eg)[i] < -geps) ) {
1150  (*ev)[i] = 0.0;
1151  }
1152  }
1153  }
1154 
1155  const ROL::Ptr<const ROL::Vector<Real> > getLowerBound(void) const {
1156  return l_;
1157  }
1158 
1159  const ROL::Ptr<const ROL::Vector<Real> > getUpperBound(void) const {
1160  return u_;
1161  }
1162 };
1163 
1164 template<class Real>
1165 class H1BoundConstraint : public ROL::BoundConstraint<Real> {
1166 private:
1167  int dim_;
1168  std::vector<Real> x_lo_;
1169  std::vector<Real> x_up_;
1170  Real min_diff_;
1171  Real scale_;
1172  ROL::Ptr<BurgersFEM<Real> > fem_;
1173  ROL::Ptr<ROL::Vector<Real> > l_;
1174  ROL::Ptr<ROL::Vector<Real> > u_;
1175 
1176  void cast_vector(ROL::Ptr<std::vector<Real> > &xvec,
1177  ROL::Vector<Real> &x) const {
1178  try {
1179  xvec = ROL::constPtrCast<std::vector<Real> >(
1180  (dynamic_cast<H1VectorPrimal<Real>&>(x)).getVector());
1181  }
1182  catch (std::exception &e) {
1183  xvec = ROL::constPtrCast<std::vector<Real> >(
1184  (dynamic_cast<H1VectorDual<Real>&>(x)).getVector());
1185  }
1186  }
1187 
1188  void cast_const_vector(ROL::Ptr<const std::vector<Real> > &xvec,
1189  const ROL::Vector<Real> &x) const {
1190  try {
1191  xvec = dynamic_cast<H1VectorPrimal<Real>&>(x).getVector();
1192  }
1193  catch (std::exception &e) {
1194  xvec = dynamic_cast<H1VectorDual<Real>&>(x).getVector();
1195  }
1196  }
1197 
1198  void axpy(std::vector<Real> &out, const Real a,
1199  const std::vector<Real> &x, const std::vector<Real> &y) const{
1200  out.resize(dim_,0.0);
1201  for (unsigned i = 0; i < dim_; i++) {
1202  out[i] = a*x[i] + y[i];
1203  }
1204  }
1205 
1206  void projection(std::vector<Real> &x) {
1207  for ( int i = 0; i < dim_; i++ ) {
1208  x[i] = std::max(x_lo_[i],std::min(x_up_[i],x[i]));
1209  }
1210  }
1211 
1212 public:
1213  H1BoundConstraint(std::vector<Real> &l, std::vector<Real> &u,
1214  const ROL::Ptr<BurgersFEM<Real> > &fem, Real scale = 1.0)
1215  : x_lo_(l), x_up_(u), scale_(scale), fem_(fem) {
1216  dim_ = x_lo_.size();
1217  for ( int i = 0; i < dim_; i++ ) {
1218  if ( i == 0 ) {
1219  min_diff_ = x_up_[i] - x_lo_[i];
1220  }
1221  else {
1222  min_diff_ = ( (min_diff_ < (x_up_[i] - x_lo_[i])) ? min_diff_ : (x_up_[i] - x_lo_[i]) );
1223  }
1224  }
1225  min_diff_ *= 0.5;
1226  l_ = ROL::makePtr<H1VectorPrimal<Real>>(
1227  ROL::makePtr<std::vector<Real>>(l), fem);
1228  u_ = ROL::makePtr<H1VectorPrimal<Real>>(
1229  ROL::makePtr<std::vector<Real>>(u), fem);
1230  }
1231 
1232  bool isFeasible( const ROL::Vector<Real> &x ) {
1233  ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1234  bool val = true;
1235  int cnt = 1;
1236  for ( int i = 0; i < dim_; i++ ) {
1237  if ( (*ex)[i] >= x_lo_[i] && (*ex)[i] <= x_up_[i] ) { cnt *= 1; }
1238  else { cnt *= 0; }
1239  }
1240  if ( cnt == 0 ) { val = false; }
1241  return val;
1242  }
1243 
1245  ROL::Ptr<std::vector<Real> > ex; cast_vector(ex,x);
1246  projection(*ex);
1247  }
1248 
1249  void pruneLowerActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &x, Real eps = Real(0)) {
1250  ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1251  ROL::Ptr<std::vector<Real> > ev; cast_vector(ev,v);
1252  Real epsn = std::min(scale_*eps,min_diff_);
1253  for ( int i = 0; i < dim_; i++ ) {
1254  if ( ((*ex)[i] <= x_lo_[i]+epsn) ) {
1255  (*ev)[i] = 0.0;
1256  }
1257  }
1258  }
1259 
1260  void pruneUpperActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &x, Real eps = Real(0)) {
1261  ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1262  ROL::Ptr<std::vector<Real> > ev; cast_vector(ev,v);
1263  Real epsn = std::min(scale_*eps,min_diff_);
1264  for ( int i = 0; i < dim_; i++ ) {
1265  if ( ((*ex)[i] >= x_up_[i]-epsn) ) {
1266  (*ev)[i] = 0.0;
1267  }
1268  }
1269  }
1270 
1271  void pruneActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &x, Real eps = Real(0)) {
1272  ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1273  ROL::Ptr<std::vector<Real> > ev; cast_vector(ev,v);
1274  Real epsn = std::min(scale_*eps,min_diff_);
1275  for ( int i = 0; i < dim_; i++ ) {
1276  if ( ((*ex)[i] <= x_lo_[i]+epsn) ||
1277  ((*ex)[i] >= x_up_[i]-epsn) ) {
1278  (*ev)[i] = 0.0;
1279  }
1280  }
1281  }
1282 
1283  void pruneLowerActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &g, const ROL::Vector<Real> &x, Real xeps = Real(0), Real geps = Real(0)) {
1284  ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1285  ROL::Ptr<const std::vector<Real> > eg; cast_const_vector(eg,g);
1286  ROL::Ptr<std::vector<Real> > ev; cast_vector(ev,v);
1287  Real epsn = std::min(scale_*xeps,min_diff_);
1288  for ( int i = 0; i < dim_; i++ ) {
1289  if ( ((*ex)[i] <= x_lo_[i]+epsn && (*eg)[i] > geps) ) {
1290  (*ev)[i] = 0.0;
1291  }
1292  }
1293  }
1294 
1295  void pruneUpperActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &g, const ROL::Vector<Real> &x, Real xeps = Real(0), Real geps = Real(0)) {
1296  ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1297  ROL::Ptr<const std::vector<Real> > eg; cast_const_vector(eg,g);
1298  ROL::Ptr<std::vector<Real> > ev; cast_vector(ev,v);
1299  Real epsn = std::min(scale_*xeps,min_diff_);
1300  for ( int i = 0; i < dim_; i++ ) {
1301  if ( ((*ex)[i] >= x_up_[i]-epsn && (*eg)[i] < -geps) ) {
1302  (*ev)[i] = 0.0;
1303  }
1304  }
1305  }
1306 
1307  void pruneActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &g, const ROL::Vector<Real> &x, Real xeps = Real(0), Real geps = Real(0)) {
1308  ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1309  ROL::Ptr<const std::vector<Real> > eg; cast_const_vector(eg,g);
1310  ROL::Ptr<std::vector<Real> > ev; cast_vector(ev,v);
1311  Real epsn = std::min(scale_*xeps,min_diff_);
1312  for ( int i = 0; i < dim_; i++ ) {
1313  if ( ((*ex)[i] <= x_lo_[i]+epsn && (*eg)[i] > geps) ||
1314  ((*ex)[i] >= x_up_[i]-epsn && (*eg)[i] < -geps) ) {
1315  (*ev)[i] = 0.0;
1316  }
1317  }
1318  }
1319 
1320  const ROL::Ptr<const ROL::Vector<Real> > getLowerBound(void) const {
1321  return l_;
1322  }
1323 
1324  const ROL::Ptr<const ROL::Vector<Real> > getUpperBound(void) const {
1325  return u_;
1326  }
1327 };
1328 
1329 template<class Real, class Ordinal>
1330 class L2VectorBatchManager : public ROL::TeuchosBatchManager<Real,Ordinal> {
1331 private:
1332  void cast_vector(ROL::Ptr<std::vector<Real> > &xvec,
1333  ROL::Vector<Real> &x) const {
1334  try {
1335  xvec = dynamic_cast<L2VectorPrimal<Real>&>(x).getVector();
1336  }
1337  catch (std::exception &e) {
1338  xvec = dynamic_cast<L2VectorDual<Real>&>(x).getVector();
1339  }
1340  }
1341 
1342 public:
1343  L2VectorBatchManager(const ROL::Ptr<const Teuchos::Comm<int> > &comm)
1344  : ROL::TeuchosBatchManager<Real,Ordinal>(comm) {}
1346  ROL::Ptr<std::vector<Real> > input_ptr;
1347  cast_vector(input_ptr,input);
1348  int dim_i = input_ptr->size();
1349  ROL::Ptr<std::vector<Real> > output_ptr;
1350  cast_vector(output_ptr,output);
1351  int dim_o = output_ptr->size();
1352  if ( dim_i != dim_o ) {
1353  std::cout << "L2VectorBatchManager: DIMENSION MISMATCH ON RANK "
1354  << ROL::TeuchosBatchManager<Real,Ordinal>::batchID() << "\n";
1355  }
1356  else {
1357  ROL::TeuchosBatchManager<Real,Ordinal>::sumAll(&(*input_ptr)[0],&(*output_ptr)[0],dim_i);
1358  }
1359  }
1360 };
1361 
1362 template<class Real, class Ordinal>
1363 class H1VectorBatchManager : public ROL::TeuchosBatchManager<Real,Ordinal> {
1364 private:
1365  void cast_vector(ROL::Ptr<std::vector<Real> > &xvec,
1366  ROL::Vector<Real> &x) const {
1367  try {
1368  xvec = dynamic_cast<H1VectorPrimal<Real>&>(x).getVector();
1369  }
1370  catch (std::exception &e) {
1371  xvec = dynamic_cast<H1VectorDual<Real>&>(x).getVector();
1372  }
1373  }
1374 
1375 public:
1376  H1VectorBatchManager(const ROL::Ptr<const Teuchos::Comm<int> > &comm)
1377  : ROL::TeuchosBatchManager<Real,Ordinal>(comm) {}
1379  ROL::Ptr<std::vector<Real> > input_ptr;
1380  cast_vector(input_ptr,input);
1381  int dim_i = input_ptr->size();
1382  ROL::Ptr<std::vector<Real> > output_ptr;
1383  cast_vector(output_ptr,output);
1384  int dim_o = output_ptr->size();
1385  if ( dim_i != dim_o ) {
1386  std::cout << "H1VectorBatchManager: DIMENSION MISMATCH ON RANK "
1387  << ROL::TeuchosBatchManager<Real,Ordinal>::batchID() << "\n";
1388  }
1389  else {
1390  ROL::TeuchosBatchManager<Real,Ordinal>::sumAll(&(*input_ptr)[0],&(*output_ptr)[0],dim_i);
1391  }
1392  }
1393 };
1394 
1395 template<class Real>
1396 Real random(const ROL::Ptr<const Teuchos::Comm<int> > &comm) {
1397  Real val = 0.0;
1398  if ( Teuchos::rank<int>(*comm)==0 ) {
1399  val = (Real)rand()/(Real)RAND_MAX;
1400  }
1401  Teuchos::broadcast<int,Real>(*comm,0,1,&val);
1402  return val;
1403 }
H1VectorPrimal< Real > DualConstraintVector
Definition: example_08.hpp:669
BurgersFEM(int nx=128, Real nl=1.0, Real cH1=1.0, Real cL2=1.0)
Definition: example_08.hpp:99
L2VectorDual(const ROL::Ptr< std::vector< Real > > &vec, const ROL::Ptr< BurgersFEM< Real > > &fem)
Definition: example_08.hpp:541
Provides the interface to evaluate simulation-based objective functions.
std::vector< Real > x_up_
Real dx_
Definition: test_04.hpp:38
ROL::Ptr< ROL::Vector< Real > > diff_
Real compute_H1_dot(const std::vector< Real > &x, const std::vector< Real > &y) const
Definition: example_08.hpp:211
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_08.hpp:748
Real cL2_
Definition: test_04.hpp:45
void apply_mass(std::vector< Real > &Mu, const std::vector< Real > &u) const
Definition: example_08.hpp:144
bool isFeasible(const ROL::Vector< Real > &x)
Check if the vector, v, is feasible.
void pruneActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &g, const ROL::Vector< Real > &x, Real xeps=Real(0), Real geps=Real(0))
void axpy(std::vector< Real > &out, const Real a, const std::vector< Real > &x, const std::vector< Real > &y) const
Definition: example_08.hpp:56
void cast_vector(ROL::Ptr< std::vector< Real > > &xvec, ROL::Vector< Real > &x) const
ROL::Ptr< L2VectorPrimal< Real > > prim_vec_
Definition: example_08.hpp:537
Real dot(const ROL::Vector< Real > &x) const
Compute where .
Definition: example_08.hpp:546
void apply_adjoint_pde_control_hessian(std::vector< Real > &ahwv, const std::vector< Real > &w, const std::vector< Real > &v, const std::vector< Real > &u, const std::vector< Real > &z)
Definition: example_08.hpp:473
ROL::Ptr< H1VectorPrimal< Real > > prim_vec_
Definition: example_08.hpp:618
ROL::Ptr< L2VectorDual< Real > > dual_vec_
Definition: test_04.hpp:503
StdVector(const Ptr< std::vector< Real >> &std_vec)
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_08.hpp:827
std::vector< Real > x_up_
void pruneUpperActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &g, const ROL::Vector< Real > &x, Real xeps=Real(0), Real geps=Real(0))
Set variables to zero if they correspond to the upper -binding set.
void apply_adjoint_pde_hessian(std::vector< Real > &ahwv, const std::vector< Real > &w, const std::vector< Real > &v, const std::vector< Real > &u, const std::vector< Real > &z) const
Definition: example_08.hpp:456
ROL::Ptr< BurgersFEM< Real > > fem_
Definition: test_04.hpp:626
H1VectorPrimal(const ROL::Ptr< std::vector< Real > > &vec, const ROL::Ptr< BurgersFEM< Real > > &fem)
Definition: example_08.hpp:585
Real u0_
Definition: test_04.hpp:41
Contains definitions of custom data types in ROL.
Ptr< const std::vector< Element > > getVector() const
Real dot(const ROL::Vector< Real > &x) const
Compute where .
Definition: example_08.hpp:627
void cast_vector(ROL::Ptr< std::vector< Real > > &xvec, ROL::Vector< Real > &x) const
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_08.hpp:802
ROL::Ptr< const std::vector< Real > > getVector() const
Definition: example_04.hpp:543
const ROL::Ptr< const ROL::Vector< Real > > getUpperBound(void) const
Return the ref count pointer to the upper bound vector.
ROL::Ptr< BurgersFEM< Real > > fem_
Definition: test_04.hpp:587
const std::vector< Real > getParameter(void) const
Real compute_H1_norm(const std::vector< Real > &r) const
Definition: example_08.hpp:231
L2BoundConstraint(std::vector< Real > &l, std::vector< Real > &u, const ROL::Ptr< BurgersFEM< Real > > &fem, Real scale=1.0)
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_08.hpp:784
void pruneLowerActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &x, Real eps=Real(0))
Set variables to zero if they correspond to the lower -active set.
H1VectorDual(const ROL::Ptr< std::vector< Real > > &vec, const ROL::Ptr< BurgersFEM< Real > > &fem)
Definition: example_08.hpp:622
bool isDualInitialized_
Definition: test_04.hpp:504
void project(ROL::Vector< Real > &x)
Project optimization variables onto the bounds.
virtual void zero()
Set to zero vector.
Definition: ROL_Vector.hpp:133
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:46
void cast_vector(ROL::Ptr< std::vector< Real > > &xvec, ROL::Vector< Real > &x) const
ROL::Ptr< BurgersFEM< Real > > fem_
H1VectorDual< Real > DualStateVector
Definition: example_08.hpp:906
H1VectorDual< Real > PrimalConstraintVector
Definition: example_08.hpp:668
void pruneLowerActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &x, Real eps=Real(0))
Set variables to zero if they correspond to the lower -active set.
void pruneLowerActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &g, const ROL::Vector< Real > &x, Real xeps=Real(0), Real geps=Real(0))
Set variables to zero if they correspond to the -binding set.
L2VectorPrimal(const ROL::Ptr< std::vector< Real > > &vec, const ROL::Ptr< BurgersFEM< Real > > &fem)
Definition: example_08.hpp:504
Real dot(const ROL::Vector< Real > &x) const
Compute where .
Definition: example_08.hpp:590
int num_dof(void) const
Definition: example_08.hpp:109
H1VectorPrimal< Real > PrimalStateVector
Definition: example_08.hpp:662
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_08.hpp:987
Constraint_BurgersControl(ROL::Ptr< BurgersFEM< Real > > &fem, bool useHessian=true)
Definition: example_08.hpp:675
Real value(const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Compute value.
Definition: example_08.hpp:925
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_08.hpp:940
Real mesh_spacing(void) const
Definition: example_08.hpp:113
ROL::Ptr< ROL::Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
Definition: example_08.hpp:560
bool isFeasible(const ROL::Vector< Real > &x)
Check if the vector, v, is feasible.
void test_inverse_mass(std::ostream &outStream=std::cout)
Definition: example_08.hpp:174
void apply_pde_jacobian(std::vector< Real > &jv, const std::vector< Real > &v, const std::vector< Real > &u, const std::vector< Real > &z) const
Definition: example_08.hpp:348
ROL::Ptr< BurgersFEM< Real > > fem_
Definition: test_04.hpp:682
const ROL::Ptr< const ROL::Vector< Real > > getLowerBound(void) const
Return the ref count pointer to the lower bound vector.
Real nl_
Definition: test_04.hpp:40
ROL::Ptr< ROL::Vector< Real > > ud_
L2VectorBatchManager(const ROL::Ptr< const Teuchos::Comm< int > > &comm)
Real cH1_
Definition: test_04.hpp:44
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_08.hpp:851
Real u1_
Definition: test_04.hpp:42
Provides the ROL::Vector interface for scalar values, to be used, for example, with scalar constraint...
void pruneUpperActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &x, Real eps=Real(0))
Set variables to zero if they correspond to the upper -active set.
ROL::Ptr< const std::vector< Real > > getVector() const
Definition: example_04.hpp:722
ROL::Ptr< ROL::Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
Definition: example_08.hpp:641
ROL::Ptr< const std::vector< Real > > getVector() const
Definition: example_04.hpp:636
ROL::Ptr< BurgersFEM< Real > > fem_
void compute_pde_jacobian(std::vector< Real > &dl, std::vector< Real > &d, std::vector< Real > &du, const std::vector< Real > &u) const
Definition: example_08.hpp:320
H1VectorDual< Real > DualStateVector
Definition: example_08.hpp:663
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_08.hpp:694
Real compute_L2_dot(const std::vector< Real > &x, const std::vector< Real > &y) const
Definition: example_08.hpp:121
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_08.hpp:982
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) const
Definition: example_08.hpp:68
std::vector< Real > x_lo_
L2VectorPrimal< Real > PrimalControlVector
Definition: example_08.hpp:665
ROL::Ptr< ROL::Vector< Real > > u_
void apply_inverse_mass(std::vector< Real > &Mu, const std::vector< Real > &u) const
Definition: example_08.hpp:161
void pruneActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &g, const ROL::Vector< Real > &x, Real xeps=Real(0), Real geps=Real(0))
void axpy(std::vector< Real > &out, const Real a, const std::vector< Real > &x, const std::vector< Real > &y) const
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_08.hpp:977
ROL::Ptr< BurgersFEM< Real > > fem_
Definition: test_04.hpp:502
ROL::Ptr< const std::vector< Real > > getVector() const
Definition: example_04.hpp:815
ROL::Ptr< BurgersFEM< Real > > fem_
Definition: test_04.hpp:541
void update(std::vector< Real > &u, const std::vector< Real > &s, const Real alpha=1.0) const
Definition: example_08.hpp:50
void pruneUpperActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &x, Real eps=Real(0))
Set variables to zero if they correspond to the upper -active set.
void projection(std::vector< Real > &x)
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_08.hpp:955
Real random(const ROL::Ptr< const Teuchos::Comm< int > > &comm)
Definition: example_05.cpp:15
void pruneActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &x, Real eps=Real(0))
int dimension() const
Return dimension of the vector space.
ROL::Ptr< ROL::Vector< Real > > l_
void apply_adjoint_control_hessian(std::vector< Real > &ahwv, const std::vector< Real > &w, const std::vector< Real > &v, const std::vector< Real > &u, const std::vector< Real > &z)
Definition: example_08.hpp:487
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_08.hpp:730
Provides the interface to apply upper and lower bound constraints.
void compute_residual(std::vector< Real > &r, const std::vector< Real > &u, const std::vector< Real > &z) const
Definition: example_08.hpp:284
void projection(std::vector< Real > &x)
void cast_const_vector(ROL::Ptr< const std::vector< Real > > &xvec, const ROL::Vector< Real > &x) const
Real nu_
Definition: test_04.hpp:39
ROL::Ptr< ROL::Vector< Real > > l_
ROL::Ptr< ROL::Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
Definition: example_08.hpp:516
H1VectorBatchManager(const ROL::Ptr< const Teuchos::Comm< int > > &comm)
void set_problem_data(const Real nu, const Real f, const Real u0, const Real u1)
Definition: example_08.hpp:102
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_08.hpp:967
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_08.hpp:712
void apply_control_jacobian(std::vector< Real > &jv, const std::vector< Real > &v, const std::vector< Real > &u, const std::vector< Real > &z) const
Definition: example_08.hpp:419
void test_inverse_H1(std::ostream &outStream=std::cout)
Definition: example_08.hpp:263
const ROL::Vector< Real > & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
Definition: example_08.hpp:602
void apply_adjoint_control_jacobian(std::vector< Real > &jv, const std::vector< Real > &v, const std::vector< Real > &u, const std::vector< Real > &z) const
Definition: example_08.hpp:430
Real dot(const ROL::Vector< Real > &x) const
Compute where .
Definition: example_08.hpp:509
void cast_vector(ROL::Ptr< std::vector< Real > > &xvec, ROL::Vector< Real > &x) const
L2VectorPrimal< Real > PrimalControlVector
Definition: example_08.hpp:908
void cast_const_vector(ROL::Ptr< const std::vector< Real > > &xvec, const ROL::Vector< Real > &x) const
H1VectorPrimal< Real > PrimalStateVector
Definition: example_08.hpp:905
std::vector< Real > x_lo_
void apply_inverse_pde_jacobian(std::vector< Real > &ijv, const std::vector< Real > &v, const std::vector< Real > &u, const std::vector< Real > &z) const
Definition: example_08.hpp:367
void axpy(std::vector< Real > &out, const Real a, const std::vector< Real > &x, const std::vector< Real > &y) const
L2VectorDual< Real > DualControlVector
Definition: example_08.hpp:666
void pruneLowerActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &g, const ROL::Vector< Real > &x, Real xeps=Real(0), Real geps=Real(0))
Set variables to zero if they correspond to the -binding set.
Real f_
Definition: test_04.hpp:43
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_08.hpp:875
void apply_H1(std::vector< Real > &Mu, const std::vector< Real > &u) const
Definition: example_08.hpp:236
const ROL::Ptr< const ROL::Vector< Real > > getLowerBound(void) const
Return the ref count pointer to the lower bound vector.
void sumAll(ROL::Vector< Real > &input, ROL::Vector< Real > &output)
void apply_inverse_H1(std::vector< Real > &Mu, const std::vector< Real > &u) const
Definition: example_08.hpp:255
void project(ROL::Vector< Real > &x)
Project optimization variables onto the bounds.
void apply_adjoint_pde_jacobian(std::vector< Real > &ajv, const std::vector< Real > &v, const std::vector< Real > &u, const std::vector< Real > &z) const
Definition: example_08.hpp:381
ROL::Ptr< BurgersFEM< Real > > fem_
ROL::Ptr< ROL::Vector< Real > > u_
ROL::Ptr< ROL::Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
Definition: example_08.hpp:597
void apply_inverse_adjoint_pde_jacobian(std::vector< Real > &iajv, const std::vector< Real > &v, const std::vector< Real > &u, const std::vector< Real > &z) const
Definition: example_08.hpp:402
void apply_adjoint_control_pde_hessian(std::vector< Real > &ahwv, const std::vector< Real > &w, const std::vector< Real > &v, const std::vector< Real > &u, const std::vector< Real > &z)
Definition: example_08.hpp:480
void pruneUpperActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &g, const ROL::Vector< Real > &x, Real xeps=Real(0), Real geps=Real(0))
Set variables to zero if they correspond to the upper -binding set.
void pruneActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &x, Real eps=Real(0))
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_08.hpp:766
int dimension() const
Return dimension of the vector space.
Definition: example_04.hpp:830
Defines the constraint operator interface for simulation-based optimization.
L2VectorDual< Real > DualControlVector
Definition: example_08.hpp:909
const ROL::Vector< Real > & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
Definition: example_08.hpp:646
bool isPrimalInitialized_
Definition: example_08.hpp:619
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_08.hpp:678
void sumAll(ROL::Vector< Real > &input, ROL::Vector< Real > &output)
bool isPrimalInitialized_
Definition: example_08.hpp:538
void scale(std::vector< Real > &u, const Real alpha=0.0) const
Definition: example_08.hpp:62
const ROL::Vector< Real > & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
Definition: example_08.hpp:565
bool isDualInitialized_
Definition: test_04.hpp:589
constexpr auto dim
ROL::Ptr< H1VectorDual< Real > > dual_vec_
Definition: test_04.hpp:588
Objective_BurgersControl(const ROL::Ptr< BurgersFEM< Real > > &fem, const ROL::Ptr< ROL::Vector< Real > > &ud, Real alpha=1.e-4)
Definition: example_08.hpp:917
const ROL::Ptr< const ROL::Vector< Real > > getUpperBound(void) const
Return the ref count pointer to the upper bound vector.
H1BoundConstraint(std::vector< Real > &l, std::vector< Real > &u, const ROL::Ptr< BurgersFEM< Real > > &fem, Real scale=1.0)
int dimension() const
Return dimension of the vector space.
Definition: example_04.hpp:651
const ROL::Vector< Real > & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
Definition: example_08.hpp:521
Real compute_L2_norm(const std::vector< Real > &r) const
Definition: example_08.hpp:139