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