ROL
example_07.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 #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  Real get_viscosity(void) const {
110  return nu_;
111  }
112 
113  int num_dof(void) const {
114  return nx_;
115  }
116 
117  Real mesh_spacing(void) const {
118  return dx_;
119  }
120 
121  /***************************************************************************/
122  /*********************** L2 INNER PRODUCT FUNCTIONS ************************/
123  /***************************************************************************/
124  // Compute L2 inner product
125  Real compute_L2_dot(const std::vector<Real> &x, const std::vector<Real> &y) const {
126  Real ip = 0.0;
127  Real c = (((int)x.size()==nx_) ? 4.0 : 2.0);
128  for (unsigned i=0; i<x.size(); i++) {
129  if ( i == 0 ) {
130  ip += dx_/6.0*(c*x[i] + x[i+1])*y[i];
131  }
132  else if ( i == x.size()-1 ) {
133  ip += dx_/6.0*(x[i-1] + c*x[i])*y[i];
134  }
135  else {
136  ip += dx_/6.0*(x[i-1] + 4.0*x[i] + x[i+1])*y[i];
137  }
138  }
139  return ip;
140  }
141 
142  // compute L2 norm
143  Real compute_L2_norm(const std::vector<Real> &r) const {
144  return std::sqrt(compute_L2_dot(r,r));
145  }
146 
147  // Apply L2 Reisz operator
148  void apply_mass(std::vector<Real> &Mu, const std::vector<Real> &u ) const {
149  Mu.resize(u.size(),0.0);
150  Real c = (((int)u.size()==nx_) ? 4.0 : 2.0);
151  for (unsigned i=0; i<u.size(); i++) {
152  if ( i == 0 ) {
153  Mu[i] = dx_/6.0*(c*u[i] + u[i+1]);
154  }
155  else if ( i == u.size()-1 ) {
156  Mu[i] = dx_/6.0*(u[i-1] + c*u[i]);
157  }
158  else {
159  Mu[i] = dx_/6.0*(u[i-1] + 4.0*u[i] + u[i+1]);
160  }
161  }
162  }
163 
164  // Apply L2 inverse Reisz operator
165  void apply_inverse_mass(std::vector<Real> &Mu, const std::vector<Real> &u) const {
166  unsigned nx = u.size();
167  // Build mass matrix
168  std::vector<Real> dl(nx-1,dx_/6.0);
169  std::vector<Real> d(nx,2.0*dx_/3.0);
170  std::vector<Real> du(nx-1,dx_/6.0);
171  if ( (int)nx != nx_ ) {
172  d[ 0] = dx_/3.0;
173  d[nx-1] = dx_/3.0;
174  }
175  linear_solve(Mu,dl,d,du,u);
176  }
177 
178  void test_inverse_mass(std::ostream &outStream = std::cout) {
179  // State Mass Matrix
180  std::vector<Real> u(nx_,0.0), Mu(nx_,0.0), iMMu(nx_,0.0), diff(nx_,0.0);
181  for (int i = 0; i < nx_; i++) {
182  u[i] = 2.0*(Real)rand()/(Real)RAND_MAX - 1.0;
183  }
184  apply_mass(Mu,u);
185  apply_inverse_mass(iMMu,Mu);
186  axpy(diff,-1.0,iMMu,u);
187  Real error = compute_L2_norm(diff);
188  Real normu = compute_L2_norm(u);
189  outStream << "Test Inverse State Mass Matrix\n";
190  outStream << " ||u - inv(M)Mu|| = " << error << "\n";
191  outStream << " ||u|| = " << normu << "\n";
192  outStream << " Relative Error = " << error/normu << "\n";
193  outStream << "\n";
194  // Control Mass Matrix
195  u.resize(nx_+2,0.0); Mu.resize(nx_+2,0.0); iMMu.resize(nx_+2,0.0); diff.resize(nx_+2,0.0);
196  for (int i = 0; i < nx_+2; i++) {
197  u[i] = 2.0*(Real)rand()/(Real)RAND_MAX - 1.0;
198  }
199  apply_mass(Mu,u);
200  apply_inverse_mass(iMMu,Mu);
201  axpy(diff,-1.0,iMMu,u);
202  error = compute_L2_norm(diff);
203  normu = compute_L2_norm(u);
204  outStream << "Test Inverse Control Mass Matrix\n";
205  outStream << " ||z - inv(M)Mz|| = " << error << "\n";
206  outStream << " ||z|| = " << normu << "\n";
207  outStream << " Relative Error = " << error/normu << "\n";
208  outStream << "\n";
209  }
210 
211  /***************************************************************************/
212  /*********************** H1 INNER PRODUCT FUNCTIONS ************************/
213  /***************************************************************************/
214  // Compute H1 inner product
215  Real compute_H1_dot(const std::vector<Real> &x, const std::vector<Real> &y) const {
216  Real ip = 0.0;
217  for (int i=0; i<nx_; i++) {
218  if ( i == 0 ) {
219  ip += cL2_*dx_/6.0*(4.0*x[i] + x[i+1])*y[i]; // Mass term
220  ip += cH1_*(2.0*x[i] - x[i+1])/dx_*y[i]; // Stiffness term
221  }
222  else if ( i == nx_-1 ) {
223  ip += cL2_*dx_/6.0*(x[i-1] + 4.0*x[i])*y[i]; // Mass term
224  ip += cH1_*(2.0*x[i] - x[i-1])/dx_*y[i]; // Stiffness term
225  }
226  else {
227  ip += cL2_*dx_/6.0*(x[i-1] + 4.0*x[i] + x[i+1])*y[i]; // Mass term
228  ip += cH1_*(2.0*x[i] - x[i-1] - x[i+1])/dx_*y[i]; // Stiffness term
229  }
230  }
231  return ip;
232  }
233 
234  // compute H1 norm
235  Real compute_H1_norm(const std::vector<Real> &r) const {
236  return std::sqrt(compute_H1_dot(r,r));
237  }
238 
239  // Apply H2 Reisz operator
240  void apply_H1(std::vector<Real> &Mu, const std::vector<Real> &u ) const {
241  Mu.resize(nx_,0.0);
242  for (int i=0; i<nx_; i++) {
243  if ( i == 0 ) {
244  Mu[i] = cL2_*dx_/6.0*(4.0*u[i] + u[i+1])
245  + cH1_*(2.0*u[i] - u[i+1])/dx_;
246  }
247  else if ( i == nx_-1 ) {
248  Mu[i] = cL2_*dx_/6.0*(u[i-1] + 4.0*u[i])
249  + cH1_*(2.0*u[i] - u[i-1])/dx_;
250  }
251  else {
252  Mu[i] = cL2_*dx_/6.0*(u[i-1] + 4.0*u[i] + u[i+1])
253  + cH1_*(2.0*u[i] - u[i-1] - u[i+1])/dx_;
254  }
255  }
256  }
257 
258  // Apply H1 inverse Reisz operator
259  void apply_inverse_H1(std::vector<Real> &Mu, const std::vector<Real> &u) const {
260  // Build mass matrix
261  std::vector<Real> dl(nx_-1,cL2_*dx_/6.0 - cH1_/dx_);
262  std::vector<Real> d(nx_,2.0*(cL2_*dx_/3.0 + cH1_/dx_));
263  std::vector<Real> du(nx_-1,cL2_*dx_/6.0 - cH1_/dx_);
264  linear_solve(Mu,dl,d,du,u);
265  }
266 
267  void test_inverse_H1(std::ostream &outStream = std::cout) {
268  std::vector<Real> u(nx_,0.0), Mu(nx_,0.0), iMMu(nx_,0.0), diff(nx_,0.0);
269  for (int i = 0; i < nx_; i++) {
270  u[i] = 2.0*(Real)rand()/(Real)RAND_MAX - 1.0;
271  }
272  apply_H1(Mu,u);
273  apply_inverse_H1(iMMu,Mu);
274  axpy(diff,-1.0,iMMu,u);
275  Real error = compute_H1_norm(diff);
276  Real normu = compute_H1_norm(u);
277  outStream << "Test Inverse State H1 Matrix\n";
278  outStream << " ||u - inv(M)Mu|| = " << error << "\n";
279  outStream << " ||u|| = " << normu << "\n";
280  outStream << " Relative Error = " << error/normu << "\n";
281  outStream << "\n";
282  }
283 
284  /***************************************************************************/
285  /*********************** PDE RESIDUAL AND SOLVE ****************************/
286  /***************************************************************************/
287  // Compute current PDE residual
288  void compute_residual(std::vector<Real> &r, const std::vector<Real> &u,
289  const std::vector<Real> &z) const {
290  r.clear();
291  r.resize(nx_,0.0);
292  for (int i=0; i<nx_; i++) {
293  // Contribution from stiffness term
294  if ( i==0 ) {
295  r[i] = nu_/dx_*(2.0*u[i]-u[i+1]);
296  }
297  else if (i==nx_-1) {
298  r[i] = nu_/dx_*(2.0*u[i]-u[i-1]);
299  }
300  else {
301  r[i] = nu_/dx_*(2.0*u[i]-u[i-1]-u[i+1]);
302  }
303  // Contribution from nonlinear term
304  if (i<nx_-1){
305  r[i] += nl_*u[i+1]*(u[i]+u[i+1])/6.0;
306  }
307  if (i>0) {
308  r[i] -= nl_*u[i-1]*(u[i-1]+u[i])/6.0;
309  }
310  // Contribution from control
311  r[i] -= dx_/6.0*(z[i]+4.0*z[i+1]+z[i+2]);
312  // Contribution from right-hand side
313  r[i] -= dx_*f_;
314  }
315  // Contribution from Dirichlet boundary terms
316  r[0] -= nl_*(u0_*u[ 0]/6.0 + u0_*u0_/6.0) + nu_*u0_/dx_;
317  r[nx_-1] += nl_*(u1_*u[nx_-1]/6.0 + u1_*u1_/6.0) - nu_*u1_/dx_;
318  }
319 
320  /***************************************************************************/
321  /*********************** PDE JACOBIAN FUNCTIONS ****************************/
322  /***************************************************************************/
323  // Build PDE Jacobian trigiagonal matrix
324  void compute_pde_jacobian(std::vector<Real> &dl, // Lower diagonal
325  std::vector<Real> &d, // Diagonal
326  std::vector<Real> &du, // Upper diagonal
327  const std::vector<Real> &u) const { // State variable
328  // Get Diagonal and Off-Diagonal Entries of linear PDE Jacobian
329  d.clear();
330  d.resize(nx_,nu_*2.0/dx_);
331  dl.clear();
332  dl.resize(nx_-1,-nu_/dx_);
333  du.clear();
334  du.resize(nx_-1,-nu_/dx_);
335  // Contribution from nonlinearity
336  for (int i=0; i<nx_; i++) {
337  if (i<nx_-1) {
338  dl[i] += nl_*(-2.0*u[i]-u[i+1])/6.0;
339  d[i] += nl_*u[i+1]/6.0;
340  }
341  if (i>0) {
342  d[i] -= nl_*u[i-1]/6.0;
343  du[i-1] += nl_*(u[i-1]+2.0*u[i])/6.0;
344  }
345  }
346  // Contribution from Dirichlet boundary conditions
347  d[ 0] -= nl_*u0_/6.0;
348  d[nx_-1] += nl_*u1_/6.0;
349  }
350 
351  // Apply PDE Jacobian to a vector
352  void apply_pde_jacobian(std::vector<Real> &jv,
353  const std::vector<Real> &v,
354  const std::vector<Real> &u,
355  const std::vector<Real> &z) const {
356  // Fill jv
357  for (int i = 0; i < nx_; i++) {
358  jv[i] = nu_/dx_*2.0*v[i];
359  if ( i > 0 ) {
360  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]);
361  }
362  if ( i < nx_-1 ) {
363  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]);
364  }
365  }
366  jv[ 0] -= nl_*u0_/6.0*v[0];
367  jv[nx_-1] += nl_*u1_/6.0*v[nx_-1];
368  }
369 
370  // Apply inverse PDE Jacobian to a vector
371  void apply_inverse_pde_jacobian(std::vector<Real> &ijv,
372  const std::vector<Real> &v,
373  const std::vector<Real> &u,
374  const std::vector<Real> &z) const {
375  // Get PDE Jacobian
376  std::vector<Real> d(nx_,0.0);
377  std::vector<Real> dl(nx_-1,0.0);
378  std::vector<Real> du(nx_-1,0.0);
379  compute_pde_jacobian(dl,d,du,u);
380  // Solve solve state sensitivity system at current time step
381  linear_solve(ijv,dl,d,du,v);
382  }
383 
384  // Apply adjoint PDE Jacobian to a vector
385  void apply_adjoint_pde_jacobian(std::vector<Real> &ajv,
386  const std::vector<Real> &v,
387  const std::vector<Real> &u,
388  const std::vector<Real> &z) const {
389  // Fill jvp
390  for (int i = 0; i < nx_; i++) {
391  ajv[i] = nu_/dx_*2.0*v[i];
392  if ( i > 0 ) {
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  if ( i < nx_-1 ) {
397  ajv[i] += -nu_/dx_*v[i+1]+nl_*(u[i+1]/6.0*v[i]
398  -(u[i+1]+2.0*u[i])/6.0*v[i+1]);
399  }
400  }
401  ajv[ 0] -= nl_*u0_/6.0*v[0];
402  ajv[nx_-1] += nl_*u1_/6.0*v[nx_-1];
403  }
404 
405  // Apply inverse adjoint PDE Jacobian to a vector
406  void apply_inverse_adjoint_pde_jacobian(std::vector<Real> &iajv,
407  const std::vector<Real> &v,
408  const std::vector<Real> &u,
409  const std::vector<Real> &z) const {
410  // Get PDE Jacobian
411  std::vector<Real> d(nx_,0.0);
412  std::vector<Real> du(nx_-1,0.0);
413  std::vector<Real> dl(nx_-1,0.0);
414  compute_pde_jacobian(dl,d,du,u);
415  // Solve solve adjoint system at current time step
416  linear_solve(iajv,dl,d,du,v,true);
417  }
418 
419  /***************************************************************************/
420  /*********************** CONTROL JACOBIAN FUNCTIONS ************************/
421  /***************************************************************************/
422  // Apply control Jacobian to a vector
423  void apply_control_jacobian(std::vector<Real> &jv,
424  const std::vector<Real> &v,
425  const std::vector<Real> &u,
426  const std::vector<Real> &z) const {
427  for (int i=0; i<nx_; i++) {
428  // Contribution from control
429  jv[i] = -dx_/6.0*(v[i]+4.0*v[i+1]+v[i+2]);
430  }
431  }
432 
433  // Apply adjoint control Jacobian to a vector
434  void apply_adjoint_control_jacobian(std::vector<Real> &jv,
435  const std::vector<Real> &v,
436  const std::vector<Real> &u,
437  const std::vector<Real> &z) const {
438  for (int i=0; i<nx_+2; i++) {
439  if ( i == 0 ) {
440  jv[i] = -dx_/6.0*v[i];
441  }
442  else if ( i == 1 ) {
443  jv[i] = -dx_/6.0*(4.0*v[i-1]+v[i]);
444  }
445  else if ( i == nx_ ) {
446  jv[i] = -dx_/6.0*(4.0*v[i-1]+v[i-2]);
447  }
448  else if ( i == nx_+1 ) {
449  jv[i] = -dx_/6.0*v[i-2];
450  }
451  else {
452  jv[i] = -dx_/6.0*(v[i-2]+4.0*v[i-1]+v[i]);
453  }
454  }
455  }
456 
457  /***************************************************************************/
458  /*********************** AJDOINT HESSIANS **********************************/
459  /***************************************************************************/
460  void apply_adjoint_pde_hessian(std::vector<Real> &ahwv,
461  const std::vector<Real> &w,
462  const std::vector<Real> &v,
463  const std::vector<Real> &u,
464  const std::vector<Real> &z) const {
465  for (int i=0; i<nx_; i++) {
466  // Contribution from nonlinear term
467  ahwv[i] = 0.0;
468  if (i<nx_-1){
469  ahwv[i] += (w[i]*v[i+1] - w[i+1]*(2.0*v[i]+v[i+1]))/6.0;
470  }
471  if (i>0) {
472  ahwv[i] += (w[i-1]*(v[i-1]+2.0*v[i]) - w[i]*v[i-1])/6.0;
473  }
474  }
475  //ahwv.assign(u.size(),0.0);
476  }
477  void apply_adjoint_pde_control_hessian(std::vector<Real> &ahwv,
478  const std::vector<Real> &w,
479  const std::vector<Real> &v,
480  const std::vector<Real> &u,
481  const std::vector<Real> &z) {
482  ahwv.assign(u.size(),0.0);
483  }
484  void apply_adjoint_control_pde_hessian(std::vector<Real> &ahwv,
485  const std::vector<Real> &w,
486  const std::vector<Real> &v,
487  const std::vector<Real> &u,
488  const std::vector<Real> &z) {
489  ahwv.assign(z.size(),0.0);
490  }
491  void apply_adjoint_control_hessian(std::vector<Real> &ahwv,
492  const std::vector<Real> &w,
493  const std::vector<Real> &v,
494  const std::vector<Real> &u,
495  const std::vector<Real> &z) {
496  ahwv.assign(z.size(),0.0);
497  }
498 };
499 
500 template<class Real>
501 class L2VectorPrimal : public ROL::Vector<Real> {
502 private:
503  ROL::Ptr<std::vector<Real> > vec_;
504  ROL::Ptr<BurgersFEM<Real> > fem_;
505 
506  mutable ROL::Ptr<L2VectorDual<Real> > dual_vec_;
507 
508 public:
509  L2VectorPrimal(const ROL::Ptr<std::vector<Real> > & vec,
510  const ROL::Ptr<BurgersFEM<Real> > &fem)
511  : vec_(vec), fem_(fem), dual_vec_(ROL::nullPtr) {}
512 
513  void set( const ROL::Vector<Real> &x ) {
514  const L2VectorPrimal &ex = dynamic_cast<const L2VectorPrimal&>(x);
515  const std::vector<Real>& xval = *ex.getVector();
516  std::copy(xval.begin(),xval.end(),vec_->begin());
517  }
518 
519  void plus( const ROL::Vector<Real> &x ) {
520  const L2VectorPrimal &ex = dynamic_cast<const L2VectorPrimal&>(x);
521  const std::vector<Real>& xval = *ex.getVector();
522  unsigned dimension = vec_->size();
523  for (unsigned i=0; i<dimension; i++) {
524  (*vec_)[i] += xval[i];
525  }
526  }
527 
528  void scale( const Real alpha ) {
529  unsigned dimension = vec_->size();
530  for (unsigned i=0; i<dimension; i++) {
531  (*vec_)[i] *= alpha;
532  }
533  }
534 
535  Real dot( const ROL::Vector<Real> &x ) const {
536  const L2VectorPrimal & ex = dynamic_cast<const L2VectorPrimal&>(x);
537  const std::vector<Real>& xval = *ex.getVector();
538  return fem_->compute_L2_dot(xval,*vec_);
539  }
540 
541  Real norm() const {
542  Real val = 0;
543  val = std::sqrt( dot(*this) );
544  return val;
545  }
546 
547  ROL::Ptr<ROL::Vector<Real> > clone() const {
548  return ROL::makePtr<L2VectorPrimal>( ROL::makePtr<std::vector<Real>>(vec_->size(),0.0),fem_);
549  }
550 
551  ROL::Ptr<const std::vector<Real> > getVector() const {
552  return vec_;
553  }
554 
555  ROL::Ptr<std::vector<Real> > getVector() {
556  return vec_;
557  }
558 
559  ROL::Ptr<ROL::Vector<Real> > basis( const int i ) const {
560  ROL::Ptr<L2VectorPrimal> e
561  = ROL::makePtr<L2VectorPrimal>( ROL::makePtr<std::vector<Real>>(vec_->size(),0.0),fem_);
562  (*e->getVector())[i] = 1.0;
563  return e;
564  }
565 
566  int dimension() const {
567  return vec_->size();
568  }
569 
570  const ROL::Vector<Real>& dual() const {
571  dual_vec_ = ROL::makePtr<L2VectorDual<Real>>(
572  ROL::makePtr<std::vector<Real>>(*vec_),fem_);
573 
574  fem_->apply_mass(*(ROL::constPtrCast<std::vector<Real> >(dual_vec_->getVector())),*vec_);
575  return *dual_vec_;
576  }
577 
578  Real apply(const ROL::Vector<Real> &x) const {
579  const L2VectorDual<Real> &ex = dynamic_cast<const L2VectorDual<Real>&>(x);
580  const std::vector<Real>& xval = *ex.getVector();
581  return std::inner_product(vec_->begin(), vec_->end(), xval.begin(), Real(0));
582  }
583 
584 };
585 
586 template<class Real>
587 class L2VectorDual : public ROL::Vector<Real> {
588 private:
589  ROL::Ptr<std::vector<Real> > vec_;
590  ROL::Ptr<BurgersFEM<Real> > fem_;
591 
592  mutable ROL::Ptr<L2VectorPrimal<Real> > dual_vec_;
593 
594 public:
595  L2VectorDual(const ROL::Ptr<std::vector<Real> > & vec,
596  const ROL::Ptr<BurgersFEM<Real> > &fem)
597  : vec_(vec), fem_(fem), dual_vec_(ROL::nullPtr) {}
598 
599  void set( const ROL::Vector<Real> &x ) {
600  const L2VectorDual &ex = dynamic_cast<const L2VectorDual&>(x);
601  const std::vector<Real>& xval = *ex.getVector();
602  std::copy(xval.begin(),xval.end(),vec_->begin());
603  }
604 
605  void plus( const ROL::Vector<Real> &x ) {
606  const L2VectorDual &ex = dynamic_cast<const L2VectorDual&>(x);
607  const std::vector<Real>& xval = *ex.getVector();
608  unsigned dimension = vec_->size();
609  for (unsigned i=0; i<dimension; i++) {
610  (*vec_)[i] += xval[i];
611  }
612  }
613 
614  void scale( const Real alpha ) {
615  unsigned dimension = vec_->size();
616  for (unsigned i=0; i<dimension; i++) {
617  (*vec_)[i] *= alpha;
618  }
619  }
620 
621  Real dot( const ROL::Vector<Real> &x ) const {
622  const L2VectorDual & ex = dynamic_cast<const L2VectorDual&>(x);
623  const std::vector<Real>& xval = *ex.getVector();
624  unsigned dimension = vec_->size();
625  std::vector<Real> Mx(dimension,0.0);
626  fem_->apply_inverse_mass(Mx,xval);
627  Real val = 0.0;
628  for (unsigned i = 0; i < dimension; i++) {
629  val += Mx[i]*(*vec_)[i];
630  }
631  return val;
632  }
633 
634  Real norm() const {
635  Real val = 0;
636  val = std::sqrt( dot(*this) );
637  return val;
638  }
639 
640  ROL::Ptr<ROL::Vector<Real> > clone() const {
641  return ROL::makePtr<L2VectorDual>( ROL::makePtr<std::vector<Real>>(vec_->size(),0.0),fem_);
642  }
643 
644  ROL::Ptr<const std::vector<Real> > getVector() const {
645  return vec_;
646  }
647 
648  ROL::Ptr<std::vector<Real> > getVector() {
649  return vec_;
650  }
651 
652  ROL::Ptr<ROL::Vector<Real> > basis( const int i ) const {
653  ROL::Ptr<L2VectorDual> e
654  = ROL::makePtr<L2VectorDual>( ROL::makePtr<std::vector<Real>>(vec_->size(),0.0),fem_);
655  (*e->getVector())[i] = 1.0;
656  return e;
657  }
658 
659  int dimension() const {
660  return vec_->size();
661  }
662 
663  const ROL::Vector<Real>& dual() const {
664  dual_vec_ = ROL::makePtr<L2VectorPrimal<Real>>(
665  ROL::makePtr<std::vector<Real>>(*vec_),fem_);
666 
667  fem_->apply_inverse_mass(*(ROL::constPtrCast<std::vector<Real> >(dual_vec_->getVector())),*vec_);
668  return *dual_vec_;
669  }
670 
671  Real apply(const ROL::Vector<Real> &x) const {
672  const L2VectorPrimal<Real> &ex = dynamic_cast<const L2VectorPrimal<Real>&>(x);
673  const std::vector<Real>& xval = *ex.getVector();
674  return std::inner_product(vec_->begin(), vec_->end(), xval.begin(), Real(0));
675  }
676 
677 };
678 
679 template<class Real>
680 class H1VectorPrimal : public ROL::Vector<Real> {
681 private:
682  ROL::Ptr<std::vector<Real> > vec_;
683  ROL::Ptr<BurgersFEM<Real> > fem_;
684 
685  mutable ROL::Ptr<H1VectorDual<Real> > dual_vec_;
686 
687 public:
688  H1VectorPrimal(const ROL::Ptr<std::vector<Real> > & vec,
689  const ROL::Ptr<BurgersFEM<Real> > &fem)
690  : vec_(vec), fem_(fem), dual_vec_(ROL::nullPtr) {}
691 
692  void set( const ROL::Vector<Real> &x ) {
693  const H1VectorPrimal &ex = dynamic_cast<const H1VectorPrimal&>(x);
694  const std::vector<Real>& xval = *ex.getVector();
695  std::copy(xval.begin(),xval.end(),vec_->begin());
696  }
697 
698  void plus( const ROL::Vector<Real> &x ) {
699  const H1VectorPrimal &ex = dynamic_cast<const H1VectorPrimal&>(x);
700  const std::vector<Real>& xval = *ex.getVector();
701  unsigned dimension = vec_->size();
702  for (unsigned i=0; i<dimension; i++) {
703  (*vec_)[i] += xval[i];
704  }
705  }
706 
707  void scale( const Real alpha ) {
708  unsigned dimension = vec_->size();
709  for (unsigned i=0; i<dimension; i++) {
710  (*vec_)[i] *= alpha;
711  }
712  }
713 
714  Real dot( const ROL::Vector<Real> &x ) const {
715  const H1VectorPrimal & ex = dynamic_cast<const H1VectorPrimal&>(x);
716  const std::vector<Real>& xval = *ex.getVector();
717  return fem_->compute_H1_dot(xval,*vec_);
718  }
719 
720  Real norm() const {
721  Real val = 0;
722  val = std::sqrt( dot(*this) );
723  return val;
724  }
725 
726  ROL::Ptr<ROL::Vector<Real> > clone() const {
727  return ROL::makePtr<H1VectorPrimal>( ROL::makePtr<std::vector<Real>>(vec_->size(),0.0),fem_);
728  }
729 
730  ROL::Ptr<const std::vector<Real> > getVector() const {
731  return vec_;
732  }
733 
734  ROL::Ptr<std::vector<Real> > getVector() {
735  return vec_;
736  }
737 
738  ROL::Ptr<ROL::Vector<Real> > basis( const int i ) const {
739  ROL::Ptr<H1VectorPrimal> e
740  = ROL::makePtr<H1VectorPrimal>( ROL::makePtr<std::vector<Real>>(vec_->size(),0.0),fem_);
741  (*e->getVector())[i] = 1.0;
742  return e;
743  }
744 
745  int dimension() const {
746  return vec_->size();
747  }
748 
749  const ROL::Vector<Real>& dual() const {
750  dual_vec_ = ROL::makePtr<H1VectorDual<Real>>(
751  ROL::makePtr<std::vector<Real>>(*vec_),fem_);
752 
753  fem_->apply_H1(*(ROL::constPtrCast<std::vector<Real> >(dual_vec_->getVector())),*vec_);
754  return *dual_vec_;
755  }
756 
757  Real apply(const ROL::Vector<Real> &x) const {
758  const H1VectorDual<Real> &ex = dynamic_cast<const H1VectorDual<Real>&>(x);
759  const std::vector<Real>& xval = *ex.getVector();
760  return std::inner_product(vec_->begin(), vec_->end(), xval.begin(), Real(0));
761  }
762 
763 };
764 
765 template<class Real>
766 class H1VectorDual : public ROL::Vector<Real> {
767 private:
768  ROL::Ptr<std::vector<Real> > vec_;
769  ROL::Ptr<BurgersFEM<Real> > fem_;
770 
771  mutable ROL::Ptr<H1VectorPrimal<Real> > dual_vec_;
772 
773 public:
774  H1VectorDual(const ROL::Ptr<std::vector<Real> > & vec,
775  const ROL::Ptr<BurgersFEM<Real> > &fem)
776  : vec_(vec), fem_(fem), dual_vec_(ROL::nullPtr) {}
777 
778  void set( const ROL::Vector<Real> &x ) {
779  const H1VectorDual &ex = dynamic_cast<const H1VectorDual&>(x);
780  const std::vector<Real>& xval = *ex.getVector();
781  std::copy(xval.begin(),xval.end(),vec_->begin());
782  }
783 
784  void plus( const ROL::Vector<Real> &x ) {
785  const H1VectorDual &ex = dynamic_cast<const H1VectorDual&>(x);
786  const std::vector<Real>& xval = *ex.getVector();
787  unsigned dimension = vec_->size();
788  for (unsigned i=0; i<dimension; i++) {
789  (*vec_)[i] += xval[i];
790  }
791  }
792 
793  void scale( const Real alpha ) {
794  unsigned dimension = vec_->size();
795  for (unsigned i=0; i<dimension; i++) {
796  (*vec_)[i] *= alpha;
797  }
798  }
799 
800  Real dot( const ROL::Vector<Real> &x ) const {
801  const H1VectorDual & ex = dynamic_cast<const H1VectorDual&>(x);
802  const std::vector<Real>& xval = *ex.getVector();
803  unsigned dimension = vec_->size();
804  std::vector<Real> Mx(dimension,0.0);
805  fem_->apply_inverse_H1(Mx,xval);
806  Real val = 0.0;
807  for (unsigned i = 0; i < dimension; i++) {
808  val += Mx[i]*(*vec_)[i];
809  }
810  return val;
811  }
812 
813  Real norm() const {
814  Real val = 0;
815  val = std::sqrt( dot(*this) );
816  return val;
817  }
818 
819  ROL::Ptr<ROL::Vector<Real> > clone() const {
820  return ROL::makePtr<H1VectorDual>( ROL::makePtr<std::vector<Real>>(vec_->size(),0.0),fem_);
821  }
822 
823  ROL::Ptr<const std::vector<Real> > getVector() const {
824  return vec_;
825  }
826 
827  ROL::Ptr<std::vector<Real> > getVector() {
828  return vec_;
829  }
830 
831  ROL::Ptr<ROL::Vector<Real> > basis( const int i ) const {
832  ROL::Ptr<H1VectorDual> e
833  = ROL::makePtr<H1VectorDual>( ROL::makePtr<std::vector<Real>>(vec_->size(),0.0),fem_);
834  (*e->getVector())[i] = 1.0;
835  return e;
836  }
837 
838  int dimension() const {
839  return vec_->size();
840  }
841 
842  const ROL::Vector<Real>& dual() const {
843  dual_vec_ = ROL::makePtr<H1VectorPrimal<Real>>(
844  ROL::makePtr<std::vector<Real>>(*vec_),fem_);
845 
846  fem_->apply_inverse_H1(*(ROL::constPtrCast<std::vector<Real> >(dual_vec_->getVector())),*vec_);
847  return *dual_vec_;
848  }
849 
850  Real apply(const ROL::Vector<Real> &x) const {
851  const H1VectorPrimal<Real> &ex = dynamic_cast<const H1VectorPrimal<Real>&>(x);
852  const std::vector<Real>& xval = *ex.getVector();
853  return std::inner_product(vec_->begin(), vec_->end(), xval.begin(), Real(0));
854  }
855 };
856 
857 
858 template<class Real>
860 private:
861 
864 
867 
870 
872 
873  ROL::Ptr<BurgersFEM<Real> > fem_;
874  bool useHessian_;
875 
876 public:
877  Constraint_BurgersControl(ROL::Ptr<BurgersFEM<Real> > &fem, bool useHessian = true)
878  : fem_(fem), useHessian_(useHessian) {}
879 
881  const ROL::Vector<Real> &z, Real &tol) {
882  ROL::Ptr<std::vector<Real> > cp =
883  dynamic_cast<PrimalConstraintVector&>(c).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_->compute_residual(*cp,*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 PrimalStateVector&>(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 
907  const std::vector<Real> param
909  fem_->set_problem_data(param[0],param[1],param[2],param[3]);
910 
911  fem_->apply_pde_jacobian(*jvp,*vp,*up,*zp);
912  }
913 
915  const ROL::Vector<Real> &z, Real &tol) {
916  ROL::Ptr<std::vector<Real> > jvp =
917  dynamic_cast<PrimalConstraintVector&>(jv).getVector();
918  ROL::Ptr<const std::vector<Real> > vp =
919  dynamic_cast<const PrimalControlVector&>(v).getVector();
920  ROL::Ptr<const std::vector<Real> > up =
921  dynamic_cast<const PrimalStateVector&>(u).getVector();
922  ROL::Ptr<const std::vector<Real> > zp =
923  dynamic_cast<const PrimalControlVector&>(z).getVector();
924 
925  const std::vector<Real> param
927  fem_->set_problem_data(param[0],param[1],param[2],param[3]);
928 
929  fem_->apply_control_jacobian(*jvp,*vp,*up,*zp);
930  }
931 
933  const ROL::Vector<Real> &z, Real &tol) {
934  ROL::Ptr<std::vector<Real> > ijvp =
935  dynamic_cast<PrimalStateVector&>(ijv).getVector();
936  ROL::Ptr<const std::vector<Real> > vp =
937  dynamic_cast<const PrimalConstraintVector&>(v).getVector();
938  ROL::Ptr<const std::vector<Real> > up =
939  dynamic_cast<const PrimalStateVector&>(u).getVector();
940  ROL::Ptr<const std::vector<Real> > zp =
941  dynamic_cast<const PrimalControlVector&>(z).getVector();
942 
943  const std::vector<Real> param
945  fem_->set_problem_data(param[0],param[1],param[2],param[3]);
946 
947  fem_->apply_inverse_pde_jacobian(*ijvp,*vp,*up,*zp);
948  }
949 
951  const ROL::Vector<Real> &z, Real &tol) {
952  ROL::Ptr<std::vector<Real> > jvp =
953  dynamic_cast<DualStateVector&>(ajv).getVector();
954  ROL::Ptr<const std::vector<Real> > vp =
955  dynamic_cast<const DualConstraintVector&>(v).getVector();
956  ROL::Ptr<const std::vector<Real> > up =
957  dynamic_cast<const PrimalStateVector&>(u).getVector();
958  ROL::Ptr<const std::vector<Real> > zp =
959  dynamic_cast<const PrimalControlVector&>(z).getVector();
960 
961  const std::vector<Real> param
963  fem_->set_problem_data(param[0],param[1],param[2],param[3]);
964 
965  fem_->apply_adjoint_pde_jacobian(*jvp,*vp,*up,*zp);
966  }
967 
969  const ROL::Vector<Real> &z, Real &tol) {
970  ROL::Ptr<std::vector<Real> > jvp =
971  dynamic_cast<DualControlVector&>(jv).getVector();
972  ROL::Ptr<const std::vector<Real> > vp =
973  dynamic_cast<const DualConstraintVector&>(v).getVector();
974  ROL::Ptr<const std::vector<Real> > up =
975  dynamic_cast<const PrimalStateVector&>(u).getVector();
976  ROL::Ptr<const std::vector<Real> > zp =
977  dynamic_cast<const PrimalControlVector&>(z).getVector();
978 
979  const std::vector<Real> param
981  fem_->set_problem_data(param[0],param[1],param[2],param[3]);
982 
983  fem_->apply_adjoint_control_jacobian(*jvp,*vp,*up,*zp);
984  }
985 
987  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
988  ROL::Ptr<std::vector<Real> > iajvp =
989  dynamic_cast<DualConstraintVector&>(iajv).getVector();
990  ROL::Ptr<const std::vector<Real> > vp =
991  dynamic_cast<const DualStateVector&>(v).getVector();
992  ROL::Ptr<const std::vector<Real> > up =
993  dynamic_cast<const PrimalStateVector&>(u).getVector();
994  ROL::Ptr<const std::vector<Real> > zp =
995  dynamic_cast<const PrimalControlVector&>(z).getVector();
996 
997  const std::vector<Real> param
999  fem_->set_problem_data(param[0],param[1],param[2],param[3]);
1000 
1001  fem_->apply_inverse_adjoint_pde_jacobian(*iajvp,*vp,*up,*zp);
1002  }
1003 
1005  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
1006  if ( useHessian_ ) {
1007  ROL::Ptr<std::vector<Real> > ahwvp =
1008  dynamic_cast<DualStateVector&>(ahwv).getVector();
1009  ROL::Ptr<const std::vector<Real> > wp =
1010  dynamic_cast<const DualConstraintVector&>(w).getVector();
1011  ROL::Ptr<const std::vector<Real> > vp =
1012  dynamic_cast<const PrimalStateVector&>(v).getVector();
1013  ROL::Ptr<const std::vector<Real> > up =
1014  dynamic_cast<const PrimalStateVector&>(u).getVector();
1015  ROL::Ptr<const std::vector<Real> > zp =
1016  dynamic_cast<const PrimalControlVector&>(z).getVector();
1017 
1018  const std::vector<Real> param
1020  fem_->set_problem_data(param[0],param[1],param[2],param[3]);
1021 
1022  fem_->apply_adjoint_pde_hessian(*ahwvp,*wp,*vp,*up,*zp);
1023  }
1024  else {
1025  ahwv.zero();
1026  }
1027  }
1028 
1030  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
1031  if ( useHessian_ ) {
1032  ROL::Ptr<std::vector<Real> > ahwvp =
1033  dynamic_cast<DualControlVector&>(ahwv).getVector();
1034  ROL::Ptr<const std::vector<Real> > wp =
1035  dynamic_cast<const DualConstraintVector&>(w).getVector();
1036  ROL::Ptr<const std::vector<Real> > vp =
1037  dynamic_cast<const PrimalStateVector&>(v).getVector();
1038  ROL::Ptr<const std::vector<Real> > up =
1039  dynamic_cast<const PrimalStateVector&>(u).getVector();
1040  ROL::Ptr<const std::vector<Real> > zp =
1041  dynamic_cast<const PrimalControlVector&>(z).getVector();
1042 
1043  const std::vector<Real> param
1045  fem_->set_problem_data(param[0],param[1],param[2],param[3]);
1046 
1047  fem_->apply_adjoint_control_pde_hessian(*ahwvp,*wp,*vp,*up,*zp);
1048  }
1049  else {
1050  ahwv.zero();
1051  }
1052  }
1054  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
1055  if ( useHessian_ ) {
1056  ROL::Ptr<std::vector<Real> > ahwvp =
1057  dynamic_cast<DualStateVector&>(ahwv).getVector();
1058  ROL::Ptr<const std::vector<Real> > wp =
1059  dynamic_cast<const DualConstraintVector&>(w).getVector();
1060  ROL::Ptr<const std::vector<Real> > vp =
1061  dynamic_cast<const PrimalControlVector&>(v).getVector();
1062  ROL::Ptr<const std::vector<Real> > up =
1063  dynamic_cast<const PrimalStateVector&>(u).getVector();
1064  ROL::Ptr<const std::vector<Real> > zp =
1065  dynamic_cast<const PrimalControlVector&>(z).getVector();
1066 
1067  const std::vector<Real> param
1069  fem_->set_problem_data(param[0],param[1],param[2],param[3]);
1070 
1071  fem_->apply_adjoint_pde_control_hessian(*ahwvp,*wp,*vp,*up,*zp);
1072  }
1073  else {
1074  ahwv.zero();
1075  }
1076  }
1078  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
1079  if ( useHessian_ ) {
1080  ROL::Ptr<std::vector<Real> > ahwvp =
1081  dynamic_cast<DualControlVector&>(ahwv).getVector();
1082  ROL::Ptr<const std::vector<Real> > wp =
1083  dynamic_cast<const DualConstraintVector&>(w).getVector();
1084  ROL::Ptr<const std::vector<Real> > vp =
1085  dynamic_cast<const PrimalControlVector&>(v).getVector();
1086  ROL::Ptr<const std::vector<Real> > up =
1087  dynamic_cast<const PrimalStateVector&>(u).getVector();
1088  ROL::Ptr<const std::vector<Real> > zp =
1089  dynamic_cast<const PrimalControlVector&>(z).getVector();
1090 
1091  const std::vector<Real> param
1093  fem_->set_problem_data(param[0],param[1],param[2],param[3]);
1094 
1095  fem_->apply_adjoint_control_hessian(*ahwvp,*wp,*vp,*up,*zp);
1096  }
1097  else {
1098  ahwv.zero();
1099  }
1100  }
1101 };
1102 
1103 template<class Real>
1104 class Objective_BurgersControl : public ROL::Objective_SimOpt<Real> {
1105 private:
1106 
1109 
1112 
1114 
1115  ROL::Ptr<BurgersFEM<Real> > fem_;
1116 
1117  Real x_;
1118  std::vector<int> indices_;
1119 
1120 public:
1122  Real x = 0.0) : fem_(fem), x_(x) {
1123  for (int i = 1; i < fem_->num_dof()+1; i++) {
1124  if ( (Real)i*(fem_->mesh_spacing()) >= x_ ) {
1125  indices_.push_back(i-1);
1126  }
1127  }
1128  }
1129 
1131 
1132  Real value( const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
1133  ROL::Ptr<const std::vector<Real> > up =
1134  dynamic_cast<const PrimalStateVector&>(u).getVector();
1135 
1136 // const std::vector<Real> param
1137 // = ROL::Objective_SimOpt<Real>::getParameter();
1138 // fem_->set_problem_data(param[0],param[1],param[2],param[3]);
1139 // Real nu = fem_->get_viscosity();
1140 //
1141 // return 0.5*nu*fem_->compute_H1_dot(*up,*up);
1142 
1143  Real val = 0.5*((((Real)indices_[0]+1.)*(fem_->mesh_spacing())-x_)
1144  *(x_+(2.-((Real)indices_[0]+1.))*(fem_->mesh_spacing()))/(fem_->mesh_spacing())
1145  +(fem_->mesh_spacing())) * (*up)[indices_[0]];
1146  for (uint i = 1; i < indices_.size(); i++) {
1147  val += (fem_->mesh_spacing())*(*up)[indices_[i]];
1148  }
1149  return -val;
1150  }
1151 
1152  void gradient_1( ROL::Vector<Real> &g, const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
1153  ROL::Ptr<std::vector<Real> > gp =
1154  dynamic_cast<DualStateVector&>(g).getVector();
1155  ROL::Ptr<const std::vector<Real> > up =
1156  dynamic_cast<const PrimalStateVector&>(u).getVector();
1157 
1158 // const std::vector<Real> param
1159 // = ROL::Objective_SimOpt<Real>::getParameter();
1160 // fem_->set_problem_data(param[0],param[1],param[2],param[3]);
1161 // Real nu = fem_->get_viscosity();
1162 //
1163 // fem_->apply_H1(*gp,*up);
1164 // g.scale(nu);
1165 
1166  g.zero();
1167  (*gp)[indices_[0]] = -0.5*((((Real)indices_[0]+1.)*(fem_->mesh_spacing())-x_)
1168  *(x_+(2.-((Real)indices_[0]+1.))*(fem_->mesh_spacing()))/(fem_->mesh_spacing())
1169  +(fem_->mesh_spacing()));
1170 
1171 
1172  for (uint i = 1; i < indices_.size(); i++) {
1173  (*gp)[indices_[i]] = -(fem_->mesh_spacing());
1174  }
1175  }
1176 
1177  void gradient_2( ROL::Vector<Real> &g, const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
1178  g.zero();
1179  }
1180 
1182  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
1183 // ROL::Ptr<std::vector<Real> > hvp =
1184 // ROL::constPtrCast<std::vector<Real> >((dynamic_cast<DualStateVector&>(hv)).getVector());
1185 // ROL::Ptr<const std::vector<Real> > vp =
1186 // (dynamic_cast<PrimalStateVector>(const_cast<ROL::Vector<Real> &&>(v))).getVector();
1187 //
1188 // const std::vector<Real> param
1189 // = ROL::Objective_SimOpt<Real>::getParameter();
1190 // fem_->set_problem_data(param[0],param[1],param[2],param[3]);
1191 // Real nu = fem_->get_viscosity();
1192 //
1193 // fem_->apply_H1(*hvp,*vp);
1194 // hv.scale(nu);
1195 
1196  hv.zero();
1197  }
1198 
1200  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
1201  hv.zero();
1202  }
1203 
1205  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
1206  hv.zero();
1207  }
1208 
1210  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
1211  hv.zero();
1212  }
1213 };
1214 
1215 template<class Real>
1216 class L2BoundConstraint : public ROL::BoundConstraint<Real> {
1217 private:
1218  int dim_;
1219  std::vector<Real> x_lo_;
1220  std::vector<Real> x_up_;
1221  Real min_diff_;
1222  Real scale_;
1223  ROL::Ptr<BurgersFEM<Real> > fem_;
1224  ROL::Ptr<ROL::Vector<Real> > l_;
1225  ROL::Ptr<ROL::Vector<Real> > u_;
1226 
1227  void cast_vector(ROL::Ptr<std::vector<Real> > &xvec,
1228  ROL::Vector<Real> &x) const {
1229  try {
1230  xvec = dynamic_cast<L2VectorPrimal<Real>&>(x).getVector();
1231  }
1232  catch (std::exception &e) {
1233  xvec = dynamic_cast<L2VectorDual<Real>&>(x).getVector();
1234  }
1235  }
1236 
1237  void cast_const_vector(ROL::Ptr<const std::vector<Real> > &xvec,
1238  const ROL::Vector<Real> &x) const {
1239  try {
1240  xvec = dynamic_cast<const L2VectorPrimal<Real>&>(x).getVector();
1241  }
1242  catch (std::exception &e) {
1243  xvec = dynamic_cast<const L2VectorDual<Real>&>(x).getVector();
1244  }
1245  }
1246 
1247  void axpy(std::vector<Real> &out, const Real a,
1248  const std::vector<Real> &x, const std::vector<Real> &y) const{
1249  out.resize(dim_,0.0);
1250  for (unsigned i = 0; i < dim_; i++) {
1251  out[i] = a*x[i] + y[i];
1252  }
1253  }
1254 
1255  void projection(std::vector<Real> &x) {
1256  for ( int i = 0; i < dim_; i++ ) {
1257  x[i] = std::max(x_lo_[i],std::min(x_up_[i],x[i]));
1258  }
1259  }
1260 
1261 public:
1262  L2BoundConstraint(std::vector<Real> &l, std::vector<Real> &u,
1263  const ROL::Ptr<BurgersFEM<Real> > &fem, Real scale = 1.0)
1264  : x_lo_(l), x_up_(u), scale_(scale), fem_(fem) {
1265  dim_ = x_lo_.size();
1266  for ( int i = 0; i < dim_; i++ ) {
1267  if ( i == 0 ) {
1268  min_diff_ = x_up_[i] - x_lo_[i];
1269  }
1270  else {
1271  min_diff_ = ( (min_diff_ < (x_up_[i] - x_lo_[i])) ? min_diff_ : (x_up_[i] - x_lo_[i]) );
1272  }
1273  }
1274  min_diff_ *= 0.5;
1275  l_ = ROL::makePtr<L2VectorPrimal<Real>>(
1276  ROL::makePtr<std::vector<Real>>(l), fem);
1277  u_ = ROL::makePtr<L2VectorPrimal<Real>>(
1278  ROL::makePtr<std::vector<Real>>(u), fem);
1279  }
1280 
1281  bool isFeasible( const ROL::Vector<Real> &x ) {
1282  ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1283  bool val = true;
1284  int cnt = 1;
1285  for ( int i = 0; i < dim_; i++ ) {
1286  if ( (*ex)[i] >= x_lo_[i] && (*ex)[i] <= x_up_[i] ) { cnt *= 1; }
1287  else { cnt *= 0; }
1288  }
1289  if ( cnt == 0 ) { val = false; }
1290  return val;
1291  }
1292 
1294  ROL::Ptr<std::vector<Real> > ex; cast_vector(ex,x);
1295  projection(*ex);
1296  }
1297 
1298  void pruneLowerActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &x, Real eps = Real(0)) {
1299  ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1300  ROL::Ptr<std::vector<Real> > ev; cast_vector(ev,v);
1301  Real epsn = std::min(scale_*eps,min_diff_);
1302  for ( int i = 0; i < dim_; i++ ) {
1303  if ( ((*ex)[i] <= x_lo_[i]+epsn) ) {
1304  (*ev)[i] = 0.0;
1305  }
1306  }
1307  }
1308 
1309  void pruneUpperActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &x, Real eps = Real(0)) {
1310  ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1311  ROL::Ptr<std::vector<Real> > ev; cast_vector(ev,v);
1312  Real epsn = std::min(scale_*eps,min_diff_);
1313  for ( int i = 0; i < dim_; i++ ) {
1314  if ( ((*ex)[i] >= x_up_[i]-epsn) ) {
1315  (*ev)[i] = 0.0;
1316  }
1317  }
1318  }
1319 
1320  void pruneActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &x, Real eps = Real(0)) {
1321  ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1322  ROL::Ptr<std::vector<Real> > ev; cast_vector(ev,v);
1323  Real epsn = std::min(scale_*eps,min_diff_);
1324  for ( int i = 0; i < dim_; i++ ) {
1325  if ( ((*ex)[i] <= x_lo_[i]+epsn) ||
1326  ((*ex)[i] >= x_up_[i]-epsn) ) {
1327  (*ev)[i] = 0.0;
1328  }
1329  }
1330  }
1331 
1332  void pruneLowerActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &g, const ROL::Vector<Real> &x, Real xeps = Real(0), Real geps = Real(0)) {
1333  ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1334  ROL::Ptr<const std::vector<Real> > eg; cast_const_vector(eg,g);
1335  ROL::Ptr<std::vector<Real> > ev; cast_vector(ev,v);
1336  Real epsn = std::min(scale_*xeps,min_diff_);
1337  for ( int i = 0; i < dim_; i++ ) {
1338  if ( ((*ex)[i] <= x_lo_[i]+epsn && (*eg)[i] > geps) ) {
1339  (*ev)[i] = 0.0;
1340  }
1341  }
1342  }
1343 
1344  void pruneUpperActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &g, const ROL::Vector<Real> &x, Real xeps = Real(0), Real geps = Real(0)) {
1345  ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1346  ROL::Ptr<const std::vector<Real> > eg; cast_const_vector(eg,g);
1347  ROL::Ptr<std::vector<Real> > ev; cast_vector(ev,v);
1348  Real epsn = std::min(scale_*xeps,min_diff_);
1349  for ( int i = 0; i < dim_; i++ ) {
1350  if ( ((*ex)[i] >= x_up_[i]-epsn && (*eg)[i] < -geps) ) {
1351  (*ev)[i] = 0.0;
1352  }
1353  }
1354  }
1355 
1356  void pruneActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &g, const ROL::Vector<Real> &x, Real xeps = Real(0), Real geps = Real(0)) {
1357  ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1358  ROL::Ptr<const std::vector<Real> > eg; cast_const_vector(eg,g);
1359  ROL::Ptr<std::vector<Real> > ev; cast_vector(ev,v);
1360  Real epsn = std::min(scale_*xeps,min_diff_);
1361  for ( int i = 0; i < dim_; i++ ) {
1362  if ( ((*ex)[i] <= x_lo_[i]+epsn && (*eg)[i] > geps) ||
1363  ((*ex)[i] >= x_up_[i]-epsn && (*eg)[i] < -geps) ) {
1364  (*ev)[i] = 0.0;
1365  }
1366  }
1367  }
1368 
1369  const ROL::Ptr<const ROL::Vector<Real> > getLowerBound(void) const {
1370  return l_;
1371  }
1372 
1373  const ROL::Ptr<const ROL::Vector<Real> > getUpperBound(void) const {
1374  return u_;
1375  }
1376 };
1377 
1378 template<class Real>
1379 class H1BoundConstraint : public ROL::BoundConstraint<Real> {
1380 private:
1381  int dim_;
1382  std::vector<Real> x_lo_;
1383  std::vector<Real> x_up_;
1384  Real min_diff_;
1385  Real scale_;
1386  ROL::Ptr<BurgersFEM<Real> > fem_;
1387  ROL::Ptr<ROL::Vector<Real> > l_;
1388  ROL::Ptr<ROL::Vector<Real> > u_;
1389 
1390  void cast_vector(ROL::Ptr<std::vector<Real> > &xvec,
1391  ROL::Vector<Real> &x) const {
1392  try {
1393  xvec = ROL::constPtrCast<std::vector<Real> >(
1394  (dynamic_cast<H1VectorPrimal<Real>&>(x)).getVector());
1395  }
1396  catch (std::exception &e) {
1397  xvec = ROL::constPtrCast<std::vector<Real> >(
1398  (dynamic_cast<H1VectorDual<Real>&>(x)).getVector());
1399  }
1400  }
1401 
1402  void cast_const_vector(ROL::Ptr<const std::vector<Real> > &xvec,
1403  const ROL::Vector<Real> &x) const {
1404  try {
1405  xvec = (dynamic_cast<H1VectorPrimal<Real>&>(
1406  const_cast<ROL::Vector<Real> &>(x))).getVector();
1407  }
1408  catch (std::exception &e) {
1409  xvec = (dynamic_cast<H1VectorDual<Real>&>(
1410  const_cast<ROL::Vector<Real> &>(x))).getVector();
1411  }
1412  }
1413 
1414  void axpy(std::vector<Real> &out, const Real a,
1415  const std::vector<Real> &x, const std::vector<Real> &y) const{
1416  out.resize(dim_,0.0);
1417  for (unsigned i = 0; i < dim_; i++) {
1418  out[i] = a*x[i] + y[i];
1419  }
1420  }
1421 
1422  void projection(std::vector<Real> &x) {
1423  for ( int i = 0; i < dim_; i++ ) {
1424  x[i] = std::max(x_lo_[i],std::min(x_up_[i],x[i]));
1425  }
1426  }
1427 
1428 public:
1429  H1BoundConstraint(std::vector<Real> &l, std::vector<Real> &u,
1430  const ROL::Ptr<BurgersFEM<Real> > &fem, Real scale = 1.0)
1431  : x_lo_(l), x_up_(u), scale_(scale), fem_(fem) {
1432  dim_ = x_lo_.size();
1433  for ( int i = 0; i < dim_; i++ ) {
1434  if ( i == 0 ) {
1435  min_diff_ = x_up_[i] - x_lo_[i];
1436  }
1437  else {
1438  min_diff_ = ( (min_diff_ < (x_up_[i] - x_lo_[i])) ? min_diff_ : (x_up_[i] - x_lo_[i]) );
1439  }
1440  }
1441  min_diff_ *= 0.5;
1442  l_ = ROL::makePtr<H1VectorPrimal<Real>>(
1443  ROL::makePtr<std::vector<Real>>(l), fem);
1444  u_ = ROL::makePtr<H1VectorPrimal<Real>>(
1445  ROL::makePtr<std::vector<Real>>(u), fem);
1446  }
1447 
1448  bool isFeasible( const ROL::Vector<Real> &x ) {
1449  ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1450  bool val = true;
1451  int cnt = 1;
1452  for ( int i = 0; i < dim_; i++ ) {
1453  if ( (*ex)[i] >= x_lo_[i] && (*ex)[i] <= x_up_[i] ) { cnt *= 1; }
1454  else { cnt *= 0; }
1455  }
1456  if ( cnt == 0 ) { val = false; }
1457  return val;
1458  }
1459 
1461  ROL::Ptr<std::vector<Real> > ex; cast_vector(ex,x);
1462  projection(*ex);
1463  }
1464 
1465  void pruneLowerActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &x, Real eps = Real(0)) {
1466  ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1467  ROL::Ptr<std::vector<Real> > ev; cast_vector(ev,v);
1468  Real epsn = std::min(scale_*eps,min_diff_);
1469  for ( int i = 0; i < dim_; i++ ) {
1470  if ( ((*ex)[i] <= x_lo_[i]+epsn) ) {
1471  (*ev)[i] = 0.0;
1472  }
1473  }
1474  }
1475 
1476  void pruneUpperActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &x, Real eps = Real(0)) {
1477  ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1478  ROL::Ptr<std::vector<Real> > ev; cast_vector(ev,v);
1479  Real epsn = std::min(scale_*eps,min_diff_);
1480  for ( int i = 0; i < dim_; i++ ) {
1481  if ( ((*ex)[i] >= x_up_[i]-epsn) ) {
1482  (*ev)[i] = 0.0;
1483  }
1484  }
1485  }
1486 
1487  void pruneActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &x, Real eps = Real(0)) {
1488  ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1489  ROL::Ptr<std::vector<Real> > ev; cast_vector(ev,v);
1490  Real epsn = std::min(scale_*eps,min_diff_);
1491  for ( int i = 0; i < dim_; i++ ) {
1492  if ( ((*ex)[i] <= x_lo_[i]+epsn) ||
1493  ((*ex)[i] >= x_up_[i]-epsn) ) {
1494  (*ev)[i] = 0.0;
1495  }
1496  }
1497  }
1498 
1499  void pruneLowerActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &g, const ROL::Vector<Real> &x, Real xeps = Real(0), Real geps = Real(0)) {
1500  ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1501  ROL::Ptr<const std::vector<Real> > eg; cast_const_vector(eg,g);
1502  ROL::Ptr<std::vector<Real> > ev; cast_vector(ev,v);
1503  Real epsn = std::min(scale_*xeps,min_diff_);
1504  for ( int i = 0; i < dim_; i++ ) {
1505  if ( ((*ex)[i] <= x_lo_[i]+epsn && (*eg)[i] > geps) ) {
1506  (*ev)[i] = 0.0;
1507  }
1508  }
1509  }
1510 
1511  void pruneUpperActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &g, const ROL::Vector<Real> &x, Real xeps = Real(0), Real geps = Real(0)) {
1512  ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1513  ROL::Ptr<const std::vector<Real> > eg; cast_const_vector(eg,g);
1514  ROL::Ptr<std::vector<Real> > ev; cast_vector(ev,v);
1515  Real epsn = std::min(scale_*xeps,min_diff_);
1516  for ( int i = 0; i < dim_; i++ ) {
1517  if ( ((*ex)[i] >= x_up_[i]-epsn && (*eg)[i] < -geps) ) {
1518  (*ev)[i] = 0.0;
1519  }
1520  }
1521  }
1522 
1523  void pruneActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &g, const ROL::Vector<Real> &x, Real xeps = Real(0), Real geps = Real(0)) {
1524  ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1525  ROL::Ptr<const std::vector<Real> > eg; cast_const_vector(eg,g);
1526  ROL::Ptr<std::vector<Real> > ev; cast_vector(ev,v);
1527  Real epsn = std::min(scale_*xeps,min_diff_);
1528  for ( int i = 0; i < dim_; i++ ) {
1529  if ( ((*ex)[i] <= x_lo_[i]+epsn && (*eg)[i] > geps) ||
1530  ((*ex)[i] >= x_up_[i]-epsn && (*eg)[i] < -geps) ) {
1531  (*ev)[i] = 0.0;
1532  }
1533  }
1534  }
1535 
1536  const ROL::Ptr<const ROL::Vector<Real> > getLowerBound(void) const {
1537  return l_;
1538  }
1539 
1540  const ROL::Ptr<const ROL::Vector<Real> > getUpperBound(void) const {
1541  return u_;
1542  }
1543 };
1544 
1545 template<class Real, class Ordinal>
1546 class L2VectorBatchManager : public ROL::TeuchosBatchManager<Real,Ordinal> {
1547 private:
1548  void cast_vector(ROL::Ptr<std::vector<Real> > &xvec,
1549  ROL::Vector<Real> &x) const {
1550  try {
1551  xvec = ROL::constPtrCast<std::vector<Real> >(
1552  (dynamic_cast<L2VectorPrimal<Real>&>(x)).getVector());
1553  }
1554  catch (std::exception &e) {
1555  xvec = ROL::constPtrCast<std::vector<Real> >(
1556  (dynamic_cast<L2VectorDual<Real>&>(x)).getVector());
1557  }
1558  }
1559 
1560 public:
1561  L2VectorBatchManager(const ROL::Ptr<const Teuchos::Comm<int> > &comm)
1562  : ROL::TeuchosBatchManager<Real,Ordinal>(comm) {}
1564  ROL::Ptr<std::vector<Real> > input_ptr;
1565  cast_vector(input_ptr,input);
1566  int dim_i = input_ptr->size();
1567  ROL::Ptr<std::vector<Real> > output_ptr;
1568  cast_vector(output_ptr,output);
1569  int dim_o = output_ptr->size();
1570  if ( dim_i != dim_o ) {
1571  std::cout << "L2VectorBatchManager: DIMENSION MISMATCH ON RANK "
1572  << ROL::TeuchosBatchManager<Real,Ordinal>::batchID() << "\n";
1573  }
1574  else {
1575  ROL::TeuchosBatchManager<Real,Ordinal>::sumAll(&(*input_ptr)[0],&(*output_ptr)[0],dim_i);
1576  }
1577  }
1578 };
1579 
1580 template<class Real, class Ordinal>
1581 class H1VectorBatchManager : public ROL::TeuchosBatchManager<Real,Ordinal> {
1582 private:
1583  void cast_vector(ROL::Ptr<std::vector<Real> > &xvec,
1584  ROL::Vector<Real> &x) const {
1585  try {
1586  xvec = ROL::constPtrCast<std::vector<Real> >(
1587  (dynamic_cast<H1VectorPrimal<Real>&>(x)).getVector());
1588  }
1589  catch (std::exception &e) {
1590  xvec = ROL::constPtrCast<std::vector<Real> >(
1591  (dynamic_cast<H1VectorDual<Real>&>(x)).getVector());
1592  }
1593  }
1594 
1595 public:
1596  H1VectorBatchManager(const ROL::Ptr<const Teuchos::Comm<int> > &comm)
1597  : ROL::TeuchosBatchManager<Real,Ordinal>(comm) {}
1599  ROL::Ptr<std::vector<Real> > input_ptr;
1600  cast_vector(input_ptr,input);
1601  int dim_i = input_ptr->size();
1602  ROL::Ptr<std::vector<Real> > output_ptr;
1603  cast_vector(output_ptr,output);
1604  int dim_o = output_ptr->size();
1605  if ( dim_i != dim_o ) {
1606  std::cout << "H1VectorBatchManager: DIMENSION MISMATCH ON RANK "
1607  << ROL::TeuchosBatchManager<Real,Ordinal>::batchID() << "\n";
1608  }
1609  else {
1610  ROL::TeuchosBatchManager<Real,Ordinal>::sumAll(&(*input_ptr)[0],&(*output_ptr)[0],dim_i);
1611  }
1612  }
1613 };
1614 
1615 template<class Real>
1616 Real random(const ROL::Ptr<const Teuchos::Comm<int> > &comm) {
1617  Real val = 0.0;
1618  if ( Teuchos::rank<int>(*comm)==0 ) {
1619  val = (Real)rand()/(Real)RAND_MAX;
1620  }
1621  Teuchos::broadcast<int,Real>(*comm,0,1,&val);
1622  return val;
1623 }
H1VectorPrimal< Real > DualConstraintVector
Definition: example_07.hpp:869
BurgersFEM(int nx=128, Real nl=1.0, Real cH1=1.0, Real cL2=1.0)
Definition: example_07.hpp:99
L2VectorDual(const ROL::Ptr< std::vector< Real > > &vec, const ROL::Ptr< BurgersFEM< Real > > &fem)
Definition: example_07.hpp:595
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
typename PV< Real >::size_type size_type
Real compute_H1_dot(const std::vector< Real > &x, const std::vector< Real > &y) const
Definition: example_07.hpp:215
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_07.hpp:950
Real norm() const
Returns where .
Definition: example_07.hpp:720
Real cL2_
Definition: test_04.hpp:45
void apply_mass(std::vector< Real > &Mu, const std::vector< Real > &u) const
Definition: example_07.hpp:148
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_07.hpp:56
ROL::Ptr< ROL::Vector< Real > > basis(const int i) const
Return i-th basis vector.
Definition: example_07.hpp:831
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_07.hpp:621
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_07.hpp:477
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
Real get_viscosity(void) const
Definition: example_07.hpp:109
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...
void plus(const ROL::Vector< Real > &x)
Compute , where .
Definition: example_07.hpp:698
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_07.hpp:460
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_07.hpp:688
Real u0_
Definition: test_04.hpp:41
Real norm() const
Returns where .
Definition: example_07.hpp:541
Contains definitions of custom data types in ROL.
Real dot(const ROL::Vector< Real > &x) const
Compute where .
Definition: example_07.hpp:800
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 ...
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_07.hpp:784
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_07.hpp:235
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_07.hpp:634
Real apply(const ROL::Vector< Real > &x) const
Apply to a dual vector. This is equivalent to the call .
Definition: example_07.hpp:671
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_07.hpp:986
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_07.hpp:774
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_07.hpp:827
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
H1VectorDual< Real > PrimalConstraintVector
Definition: example_07.hpp:868
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_07.hpp:850
void scale(const Real alpha)
Compute where .
Definition: example_07.hpp:528
std::vector< Real >::size_type uint
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_07.hpp:509
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_07.hpp:714
int num_dof(void) const
Definition: example_07.hpp:113
H1VectorPrimal< Real > PrimalStateVector
Definition: example_07.hpp:862
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_07.hpp:877
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_07.hpp:519
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_07.hpp:117
ROL::Ptr< ROL::Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
Definition: example_07.hpp:640
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_07.hpp:178
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_07.hpp:352
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
L2VectorBatchManager(const ROL::Ptr< const Teuchos::Comm< int > > &comm)
ROL::Ptr< ROL::Vector< Real > > basis(const int i) const
Return i-th basis vector.
Definition: example_07.hpp:559
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_07.hpp:819
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_07.hpp:692
void scale(const Real alpha)
Compute where .
Definition: example_07.hpp:614
void compute_pde_jacobian(std::vector< Real > &dl, std::vector< Real > &d, std::vector< Real > &du, const std::vector< Real > &u) const
Definition: example_07.hpp:324
ROL::Ptr< ROL::Vector< Real > > basis(const int i) const
Return i-th basis vector.
Definition: example_07.hpp:652
H1VectorDual< Real > DualStateVector
Definition: example_07.hpp:863
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_07.hpp:896
Real compute_L2_dot(const std::vector< Real > &x, const std::vector< Real > &y) const
Definition: example_07.hpp:125
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_07.hpp:68
std::vector< Real > x_lo_
ROL::Ptr< std::vector< Real > > getVector()
Definition: example_07.hpp:734
L2VectorPrimal< Real > PrimalControlVector
Definition: example_07.hpp:865
ROL::Ptr< ROL::Vector< Real > > u_
void apply_inverse_mass(std::vector< Real > &Mu, const std::vector< Real > &u) const
Definition: example_07.hpp:165
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_07.hpp:555
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_07.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 scale(const Real alpha)
Compute where .
Definition: example_07.hpp:707
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.
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))
void set(const ROL::Vector< Real > &x)
Set where .
Definition: example_07.hpp:778
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_07.hpp:757
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_07.hpp:491
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_07.hpp:932
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_07.hpp:288
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_07.hpp:547
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_07.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.
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_07.hpp:914
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_07.hpp:423
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_07.hpp:267
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_07.hpp:749
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_07.hpp:434
Real dot(const ROL::Vector< Real > &x) const
Compute where .
Definition: example_07.hpp:535
void cast_vector(ROL::Ptr< std::vector< Real > > &xvec, ROL::Vector< Real > &x) const
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_07.hpp:371
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_07.hpp:866
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
std::vector< int > indices_
ROL::Ptr< ROL::Vector< Real > > basis(const int i) const
Return i-th basis vector.
Definition: example_07.hpp:738
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_07.hpp:240
void plus(const ROL::Vector< Real > &x)
Compute , where .
Definition: example_07.hpp:605
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_07.hpp:259
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_07.hpp:385
void set(const ROL::Vector< Real > &x)
Set where .
Definition: example_07.hpp:599
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_07.hpp:578
ROL::Ptr< ROL::Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
Definition: example_07.hpp:726
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_07.hpp:406
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_07.hpp:484
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.
Objective_BurgersControl(const ROL::Ptr< BurgersFEM< Real > > &fem, Real x=0.0)
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_07.hpp:968
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_07.hpp:842
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_07.hpp:880
void sumAll(ROL::Vector< Real > &input, ROL::Vector< Real > &output)
void set(const ROL::Vector< Real > &x)
Set where .
Definition: example_07.hpp:513
Real norm() const
Returns where .
Definition: example_07.hpp:813
void scale(std::vector< Real > &u, const Real alpha=0.0) const
Definition: example_07.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_07.hpp:663
constexpr auto dim
ROL::Ptr< H1VectorDual< Real > > dual_vec_
Definition: test_04.hpp:588
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_07.hpp:648
std::vector< Real >::size_type uint
Definition: example_07.hpp:871
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_07.hpp:793
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_07.hpp:570
Real compute_L2_norm(const std::vector< Real > &r) const
Definition: example_07.hpp:143
ROL::Ptr< L2VectorPrimal< Real > > dual_vec_
Definition: test_04.hpp:542