ROL
example_07.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ************************************************************************
3 //
4 // Rapid Optimization Library (ROL) Package
5 // Copyright (2014) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact lead developers:
38 // Drew Kouri (dpkouri@sandia.gov) and
39 // Denis Ridzal (dridzal@sandia.gov)
40 //
41 // ************************************************************************
42 // @HEADER
43 
49 #include "ROL_Types.hpp"
50 #include "ROL_Vector.hpp"
51 #include "ROL_BoundConstraint.hpp"
53 #include "ROL_Objective_SimOpt.hpp"
55 
56 #include "Teuchos_LAPACK.hpp"
57 
58 template<class Real>
59 class L2VectorPrimal;
60 
61 template<class Real>
62 class L2VectorDual;
63 
64 template<class Real>
65 class H1VectorPrimal;
66 
67 template<class Real>
68 class H1VectorDual;
69 
70 template<class Real>
71 class BurgersFEM {
72 private:
73  int nx_;
74  Real dx_;
75  Real nu_;
76  Real nl_;
77  Real u0_;
78  Real u1_;
79  Real f_;
80  Real cH1_;
81  Real cL2_;
82 
83 private:
84  void update(std::vector<Real> &u, const std::vector<Real> &s, const Real alpha=1.0) const {
85  for (unsigned i=0; i<u.size(); i++) {
86  u[i] += alpha*s[i];
87  }
88  }
89 
90  void axpy(std::vector<Real> &out, const Real a, const std::vector<Real> &x, const std::vector<Real> &y) const {
91  for (unsigned i=0; i < x.size(); i++) {
92  out[i] = a*x[i] + y[i];
93  }
94  }
95 
96  void scale(std::vector<Real> &u, const Real alpha=0.0) const {
97  for (unsigned i=0; i<u.size(); i++) {
98  u[i] *= alpha;
99  }
100  }
101 
102  void linear_solve(std::vector<Real> &u, std::vector<Real> &dl, std::vector<Real> &d, std::vector<Real> &du,
103  const std::vector<Real> &r, const bool transpose = false) const {
104  if ( r.size() == 1 ) {
105  u.resize(1,r[0]/d[0]);
106  }
107  else if ( r.size() == 2 ) {
108  u.resize(2,0.0);
109  Real det = d[0]*d[1] - dl[0]*du[0];
110  u[0] = (d[1]*r[0] - du[0]*r[1])/det;
111  u[1] = (d[0]*r[1] - dl[0]*r[0])/det;
112  }
113  else {
114  u.assign(r.begin(),r.end());
115  // Perform LDL factorization
116  Teuchos::LAPACK<int,Real> lp;
117  std::vector<Real> du2(r.size()-2,0.0);
118  std::vector<int> ipiv(r.size(),0);
119  int info;
120  int dim = r.size();
121  int ldb = r.size();
122  int nhrs = 1;
123  lp.GTTRF(dim,&dl[0],&d[0],&du[0],&du2[0],&ipiv[0],&info);
124  char trans = 'N';
125  if ( transpose ) {
126  trans = 'T';
127  }
128  lp.GTTRS(trans,dim,nhrs,&dl[0],&d[0],&du[0],&du2[0],&ipiv[0],&u[0],ldb,&info);
129  }
130  }
131 
132 public:
133  BurgersFEM(int nx = 128, Real nl = 1.0, Real cH1 = 1.0, Real cL2 = 1.0)
134  : nx_(nx), dx_(1.0/((Real)nx+1.0)), nl_(nl), cH1_(cH1), cL2_(cL2) {}
135 
136  void set_problem_data(const Real nu, const Real f, const Real u0, const Real u1) {
137  nu_ = std::pow(10.0,nu-2.0);
138  f_ = 0.01*f;
139  u0_ = 1.0+0.001*u0;
140  u1_ = 0.001*u1;
141  }
142 
143  Real get_viscosity(void) const {
144  return nu_;
145  }
146 
147  int num_dof(void) const {
148  return nx_;
149  }
150 
151  Real mesh_spacing(void) const {
152  return dx_;
153  }
154 
155  /***************************************************************************/
156  /*********************** L2 INNER PRODUCT FUNCTIONS ************************/
157  /***************************************************************************/
158  // Compute L2 inner product
159  Real compute_L2_dot(const std::vector<Real> &x, const std::vector<Real> &y) const {
160  Real ip = 0.0;
161  Real c = (((int)x.size()==nx_) ? 4.0 : 2.0);
162  for (unsigned i=0; i<x.size(); i++) {
163  if ( i == 0 ) {
164  ip += dx_/6.0*(c*x[i] + x[i+1])*y[i];
165  }
166  else if ( i == x.size()-1 ) {
167  ip += dx_/6.0*(x[i-1] + c*x[i])*y[i];
168  }
169  else {
170  ip += dx_/6.0*(x[i-1] + 4.0*x[i] + x[i+1])*y[i];
171  }
172  }
173  return ip;
174  }
175 
176  // compute L2 norm
177  Real compute_L2_norm(const std::vector<Real> &r) const {
178  return std::sqrt(compute_L2_dot(r,r));
179  }
180 
181  // Apply L2 Reisz operator
182  void apply_mass(std::vector<Real> &Mu, const std::vector<Real> &u ) const {
183  Mu.resize(u.size(),0.0);
184  Real c = (((int)u.size()==nx_) ? 4.0 : 2.0);
185  for (unsigned i=0; i<u.size(); i++) {
186  if ( i == 0 ) {
187  Mu[i] = dx_/6.0*(c*u[i] + u[i+1]);
188  }
189  else if ( i == u.size()-1 ) {
190  Mu[i] = dx_/6.0*(u[i-1] + c*u[i]);
191  }
192  else {
193  Mu[i] = dx_/6.0*(u[i-1] + 4.0*u[i] + u[i+1]);
194  }
195  }
196  }
197 
198  // Apply L2 inverse Reisz operator
199  void apply_inverse_mass(std::vector<Real> &Mu, const std::vector<Real> &u) const {
200  unsigned nx = u.size();
201  // Build mass matrix
202  std::vector<Real> dl(nx-1,dx_/6.0);
203  std::vector<Real> d(nx,2.0*dx_/3.0);
204  std::vector<Real> du(nx-1,dx_/6.0);
205  if ( (int)nx != nx_ ) {
206  d[ 0] = dx_/3.0;
207  d[nx-1] = dx_/3.0;
208  }
209  linear_solve(Mu,dl,d,du,u);
210  }
211 
212  void test_inverse_mass(std::ostream &outStream = std::cout) {
213  // State Mass Matrix
214  std::vector<Real> u(nx_,0.0), Mu(nx_,0.0), iMMu(nx_,0.0), diff(nx_,0.0);
215  for (int i = 0; i < nx_; i++) {
216  u[i] = 2.0*(Real)rand()/(Real)RAND_MAX - 1.0;
217  }
218  apply_mass(Mu,u);
219  apply_inverse_mass(iMMu,Mu);
220  axpy(diff,-1.0,iMMu,u);
221  Real error = compute_L2_norm(diff);
222  Real normu = compute_L2_norm(u);
223  outStream << "Test Inverse State Mass Matrix\n";
224  outStream << " ||u - inv(M)Mu|| = " << error << "\n";
225  outStream << " ||u|| = " << normu << "\n";
226  outStream << " Relative Error = " << error/normu << "\n";
227  outStream << "\n";
228  // Control Mass Matrix
229  u.resize(nx_+2,0.0); Mu.resize(nx_+2,0.0); iMMu.resize(nx_+2,0.0); diff.resize(nx_+2,0.0);
230  for (int i = 0; i < nx_+2; i++) {
231  u[i] = 2.0*(Real)rand()/(Real)RAND_MAX - 1.0;
232  }
233  apply_mass(Mu,u);
234  apply_inverse_mass(iMMu,Mu);
235  axpy(diff,-1.0,iMMu,u);
236  error = compute_L2_norm(diff);
237  normu = compute_L2_norm(u);
238  outStream << "Test Inverse Control Mass Matrix\n";
239  outStream << " ||z - inv(M)Mz|| = " << error << "\n";
240  outStream << " ||z|| = " << normu << "\n";
241  outStream << " Relative Error = " << error/normu << "\n";
242  outStream << "\n";
243  }
244 
245  /***************************************************************************/
246  /*********************** H1 INNER PRODUCT FUNCTIONS ************************/
247  /***************************************************************************/
248  // Compute H1 inner product
249  Real compute_H1_dot(const std::vector<Real> &x, const std::vector<Real> &y) const {
250  Real ip = 0.0;
251  for (int i=0; i<nx_; i++) {
252  if ( i == 0 ) {
253  ip += cL2_*dx_/6.0*(4.0*x[i] + x[i+1])*y[i]; // Mass term
254  ip += cH1_*(2.0*x[i] - x[i+1])/dx_*y[i]; // Stiffness term
255  }
256  else if ( i == nx_-1 ) {
257  ip += cL2_*dx_/6.0*(x[i-1] + 4.0*x[i])*y[i]; // Mass term
258  ip += cH1_*(2.0*x[i] - x[i-1])/dx_*y[i]; // Stiffness term
259  }
260  else {
261  ip += cL2_*dx_/6.0*(x[i-1] + 4.0*x[i] + x[i+1])*y[i]; // Mass term
262  ip += cH1_*(2.0*x[i] - x[i-1] - x[i+1])/dx_*y[i]; // Stiffness term
263  }
264  }
265  return ip;
266  }
267 
268  // compute H1 norm
269  Real compute_H1_norm(const std::vector<Real> &r) const {
270  return std::sqrt(compute_H1_dot(r,r));
271  }
272 
273  // Apply H2 Reisz operator
274  void apply_H1(std::vector<Real> &Mu, const std::vector<Real> &u ) const {
275  Mu.resize(nx_,0.0);
276  for (int i=0; i<nx_; i++) {
277  if ( i == 0 ) {
278  Mu[i] = cL2_*dx_/6.0*(4.0*u[i] + u[i+1])
279  + cH1_*(2.0*u[i] - u[i+1])/dx_;
280  }
281  else if ( i == nx_-1 ) {
282  Mu[i] = cL2_*dx_/6.0*(u[i-1] + 4.0*u[i])
283  + cH1_*(2.0*u[i] - u[i-1])/dx_;
284  }
285  else {
286  Mu[i] = cL2_*dx_/6.0*(u[i-1] + 4.0*u[i] + u[i+1])
287  + cH1_*(2.0*u[i] - u[i-1] - u[i+1])/dx_;
288  }
289  }
290  }
291 
292  // Apply H1 inverse Reisz operator
293  void apply_inverse_H1(std::vector<Real> &Mu, const std::vector<Real> &u) const {
294  // Build mass matrix
295  std::vector<Real> dl(nx_-1,cL2_*dx_/6.0 - cH1_/dx_);
296  std::vector<Real> d(nx_,2.0*(cL2_*dx_/3.0 + cH1_/dx_));
297  std::vector<Real> du(nx_-1,cL2_*dx_/6.0 - cH1_/dx_);
298  linear_solve(Mu,dl,d,du,u);
299  }
300 
301  void test_inverse_H1(std::ostream &outStream = std::cout) {
302  std::vector<Real> u(nx_,0.0), Mu(nx_,0.0), iMMu(nx_,0.0), diff(nx_,0.0);
303  for (int i = 0; i < nx_; i++) {
304  u[i] = 2.0*(Real)rand()/(Real)RAND_MAX - 1.0;
305  }
306  apply_H1(Mu,u);
307  apply_inverse_H1(iMMu,Mu);
308  axpy(diff,-1.0,iMMu,u);
309  Real error = compute_H1_norm(diff);
310  Real normu = compute_H1_norm(u);
311  outStream << "Test Inverse State H1 Matrix\n";
312  outStream << " ||u - inv(M)Mu|| = " << error << "\n";
313  outStream << " ||u|| = " << normu << "\n";
314  outStream << " Relative Error = " << error/normu << "\n";
315  outStream << "\n";
316  }
317 
318  /***************************************************************************/
319  /*********************** PDE RESIDUAL AND SOLVE ****************************/
320  /***************************************************************************/
321  // Compute current PDE residual
322  void compute_residual(std::vector<Real> &r, const std::vector<Real> &u,
323  const std::vector<Real> &z) const {
324  r.clear();
325  r.resize(nx_,0.0);
326  for (int i=0; i<nx_; i++) {
327  // Contribution from stiffness term
328  if ( i==0 ) {
329  r[i] = nu_/dx_*(2.0*u[i]-u[i+1]);
330  }
331  else if (i==nx_-1) {
332  r[i] = nu_/dx_*(2.0*u[i]-u[i-1]);
333  }
334  else {
335  r[i] = nu_/dx_*(2.0*u[i]-u[i-1]-u[i+1]);
336  }
337  // Contribution from nonlinear term
338  if (i<nx_-1){
339  r[i] += nl_*u[i+1]*(u[i]+u[i+1])/6.0;
340  }
341  if (i>0) {
342  r[i] -= nl_*u[i-1]*(u[i-1]+u[i])/6.0;
343  }
344  // Contribution from control
345  r[i] -= dx_/6.0*(z[i]+4.0*z[i+1]+z[i+2]);
346  // Contribution from right-hand side
347  r[i] -= dx_*f_;
348  }
349  // Contribution from Dirichlet boundary terms
350  r[0] -= nl_*(u0_*u[ 0]/6.0 + u0_*u0_/6.0) + nu_*u0_/dx_;
351  r[nx_-1] += nl_*(u1_*u[nx_-1]/6.0 + u1_*u1_/6.0) - nu_*u1_/dx_;
352  }
353 
354  /***************************************************************************/
355  /*********************** PDE JACOBIAN FUNCTIONS ****************************/
356  /***************************************************************************/
357  // Build PDE Jacobian trigiagonal matrix
358  void compute_pde_jacobian(std::vector<Real> &dl, // Lower diagonal
359  std::vector<Real> &d, // Diagonal
360  std::vector<Real> &du, // Upper diagonal
361  const std::vector<Real> &u) const { // State variable
362  // Get Diagonal and Off-Diagonal Entries of linear PDE Jacobian
363  d.clear();
364  d.resize(nx_,nu_*2.0/dx_);
365  dl.clear();
366  dl.resize(nx_-1,-nu_/dx_);
367  du.clear();
368  du.resize(nx_-1,-nu_/dx_);
369  // Contribution from nonlinearity
370  for (int i=0; i<nx_; i++) {
371  if (i<nx_-1) {
372  dl[i] += nl_*(-2.0*u[i]-u[i+1])/6.0;
373  d[i] += nl_*u[i+1]/6.0;
374  }
375  if (i>0) {
376  d[i] -= nl_*u[i-1]/6.0;
377  du[i-1] += nl_*(u[i-1]+2.0*u[i])/6.0;
378  }
379  }
380  // Contribution from Dirichlet boundary conditions
381  d[ 0] -= nl_*u0_/6.0;
382  d[nx_-1] += nl_*u1_/6.0;
383  }
384 
385  // Apply PDE Jacobian to a vector
386  void apply_pde_jacobian(std::vector<Real> &jv,
387  const std::vector<Real> &v,
388  const std::vector<Real> &u,
389  const std::vector<Real> &z) const {
390  // Fill jv
391  for (int i = 0; i < nx_; i++) {
392  jv[i] = nu_/dx_*2.0*v[i];
393  if ( i > 0 ) {
394  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]);
395  }
396  if ( i < nx_-1 ) {
397  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]);
398  }
399  }
400  jv[ 0] -= nl_*u0_/6.0*v[0];
401  jv[nx_-1] += nl_*u1_/6.0*v[nx_-1];
402  }
403 
404  // Apply inverse PDE Jacobian to a vector
405  void apply_inverse_pde_jacobian(std::vector<Real> &ijv,
406  const std::vector<Real> &v,
407  const std::vector<Real> &u,
408  const std::vector<Real> &z) const {
409  // Get PDE Jacobian
410  std::vector<Real> d(nx_,0.0);
411  std::vector<Real> dl(nx_-1,0.0);
412  std::vector<Real> du(nx_-1,0.0);
413  compute_pde_jacobian(dl,d,du,u);
414  // Solve solve state sensitivity system at current time step
415  linear_solve(ijv,dl,d,du,v);
416  }
417 
418  // Apply adjoint PDE Jacobian to a vector
419  void apply_adjoint_pde_jacobian(std::vector<Real> &ajv,
420  const std::vector<Real> &v,
421  const std::vector<Real> &u,
422  const std::vector<Real> &z) const {
423  // Fill jvp
424  for (int i = 0; i < nx_; i++) {
425  ajv[i] = nu_/dx_*2.0*v[i];
426  if ( i > 0 ) {
427  ajv[i] += -nu_/dx_*v[i-1]-nl_*(u[i-1]/6.0*v[i]
428  -(u[i-1]+2.0*u[i])/6.0*v[i-1]);
429  }
430  if ( i < nx_-1 ) {
431  ajv[i] += -nu_/dx_*v[i+1]+nl_*(u[i+1]/6.0*v[i]
432  -(u[i+1]+2.0*u[i])/6.0*v[i+1]);
433  }
434  }
435  ajv[ 0] -= nl_*u0_/6.0*v[0];
436  ajv[nx_-1] += nl_*u1_/6.0*v[nx_-1];
437  }
438 
439  // Apply inverse adjoint PDE Jacobian to a vector
440  void apply_inverse_adjoint_pde_jacobian(std::vector<Real> &iajv,
441  const std::vector<Real> &v,
442  const std::vector<Real> &u,
443  const std::vector<Real> &z) const {
444  // Get PDE Jacobian
445  std::vector<Real> d(nx_,0.0);
446  std::vector<Real> du(nx_-1,0.0);
447  std::vector<Real> dl(nx_-1,0.0);
448  compute_pde_jacobian(dl,d,du,u);
449  // Solve solve adjoint system at current time step
450  linear_solve(iajv,dl,d,du,v,true);
451  }
452 
453  /***************************************************************************/
454  /*********************** CONTROL JACOBIAN FUNCTIONS ************************/
455  /***************************************************************************/
456  // Apply control Jacobian to a vector
457  void apply_control_jacobian(std::vector<Real> &jv,
458  const std::vector<Real> &v,
459  const std::vector<Real> &u,
460  const std::vector<Real> &z) const {
461  for (int i=0; i<nx_; i++) {
462  // Contribution from control
463  jv[i] = -dx_/6.0*(v[i]+4.0*v[i+1]+v[i+2]);
464  }
465  }
466 
467  // Apply adjoint control Jacobian to a vector
468  void apply_adjoint_control_jacobian(std::vector<Real> &jv,
469  const std::vector<Real> &v,
470  const std::vector<Real> &u,
471  const std::vector<Real> &z) const {
472  for (int i=0; i<nx_+2; i++) {
473  if ( i == 0 ) {
474  jv[i] = -dx_/6.0*v[i];
475  }
476  else if ( i == 1 ) {
477  jv[i] = -dx_/6.0*(4.0*v[i-1]+v[i]);
478  }
479  else if ( i == nx_ ) {
480  jv[i] = -dx_/6.0*(4.0*v[i-1]+v[i-2]);
481  }
482  else if ( i == nx_+1 ) {
483  jv[i] = -dx_/6.0*v[i-2];
484  }
485  else {
486  jv[i] = -dx_/6.0*(v[i-2]+4.0*v[i-1]+v[i]);
487  }
488  }
489  }
490 
491  /***************************************************************************/
492  /*********************** AJDOINT HESSIANS **********************************/
493  /***************************************************************************/
494  void apply_adjoint_pde_hessian(std::vector<Real> &ahwv,
495  const std::vector<Real> &w,
496  const std::vector<Real> &v,
497  const std::vector<Real> &u,
498  const std::vector<Real> &z) const {
499  for (int i=0; i<nx_; i++) {
500  // Contribution from nonlinear term
501  ahwv[i] = 0.0;
502  if (i<nx_-1){
503  ahwv[i] += (w[i]*v[i+1] - w[i+1]*(2.0*v[i]+v[i+1]))/6.0;
504  }
505  if (i>0) {
506  ahwv[i] += (w[i-1]*(v[i-1]+2.0*v[i]) - w[i]*v[i-1])/6.0;
507  }
508  }
509  //ahwv.assign(u.size(),0.0);
510  }
511  void apply_adjoint_pde_control_hessian(std::vector<Real> &ahwv,
512  const std::vector<Real> &w,
513  const std::vector<Real> &v,
514  const std::vector<Real> &u,
515  const std::vector<Real> &z) {
516  ahwv.assign(u.size(),0.0);
517  }
518  void apply_adjoint_control_pde_hessian(std::vector<Real> &ahwv,
519  const std::vector<Real> &w,
520  const std::vector<Real> &v,
521  const std::vector<Real> &u,
522  const std::vector<Real> &z) {
523  ahwv.assign(z.size(),0.0);
524  }
525  void apply_adjoint_control_hessian(std::vector<Real> &ahwv,
526  const std::vector<Real> &w,
527  const std::vector<Real> &v,
528  const std::vector<Real> &u,
529  const std::vector<Real> &z) {
530  ahwv.assign(z.size(),0.0);
531  }
532 };
533 
534 template<class Real>
535 class L2VectorPrimal : public ROL::Vector<Real> {
536 private:
537  Teuchos::RCP<std::vector<Real> > vec_;
538  Teuchos::RCP<BurgersFEM<Real> > fem_;
539 
540  mutable Teuchos::RCP<L2VectorDual<Real> > dual_vec_;
541 
542 public:
543  L2VectorPrimal(const Teuchos::RCP<std::vector<Real> > & vec,
544  const Teuchos::RCP<BurgersFEM<Real> > &fem)
545  : vec_(vec), fem_(fem), dual_vec_(Teuchos::null) {}
546 
547  void set( const ROL::Vector<Real> &x ) {
548  const L2VectorPrimal &ex = Teuchos::dyn_cast<const L2VectorPrimal>(x);
549  const std::vector<Real>& xval = *ex.getVector();
550  std::copy(xval.begin(),xval.end(),vec_->begin());
551  }
552 
553  void plus( const ROL::Vector<Real> &x ) {
554  const L2VectorPrimal &ex = Teuchos::dyn_cast<const L2VectorPrimal>(x);
555  const std::vector<Real>& xval = *ex.getVector();
556  unsigned dimension = vec_->size();
557  for (unsigned i=0; i<dimension; i++) {
558  (*vec_)[i] += xval[i];
559  }
560  }
561 
562  void scale( const Real alpha ) {
563  unsigned dimension = vec_->size();
564  for (unsigned i=0; i<dimension; i++) {
565  (*vec_)[i] *= alpha;
566  }
567  }
568 
569  Real dot( const ROL::Vector<Real> &x ) const {
570  const L2VectorPrimal & ex = Teuchos::dyn_cast<const L2VectorPrimal>(x);
571  const std::vector<Real>& xval = *ex.getVector();
572  return fem_->compute_L2_dot(xval,*vec_);
573  }
574 
575  Real norm() const {
576  Real val = 0;
577  val = std::sqrt( dot(*this) );
578  return val;
579  }
580 
581  Teuchos::RCP<ROL::Vector<Real> > clone() const {
582  return Teuchos::rcp( new L2VectorPrimal( Teuchos::rcp(new std::vector<Real>(vec_->size(),0.0)),fem_));
583  }
584 
585  Teuchos::RCP<const std::vector<Real> > getVector() const {
586  return vec_;
587  }
588 
589  Teuchos::RCP<std::vector<Real> > getVector() {
590  return vec_;
591  }
592 
593  Teuchos::RCP<ROL::Vector<Real> > basis( const int i ) const {
594  Teuchos::RCP<L2VectorPrimal> e
595  = Teuchos::rcp( new L2VectorPrimal( Teuchos::rcp(new std::vector<Real>(vec_->size(),0.0)),fem_));
596  (*e->getVector())[i] = 1.0;
597  return e;
598  }
599 
600  int dimension() const {
601  return vec_->size();
602  }
603 
604  const ROL::Vector<Real>& dual() const {
605  dual_vec_ = Teuchos::rcp(new L2VectorDual<Real>(
606  Teuchos::rcp(new std::vector<Real>(*vec_)),fem_));
607 
608  fem_->apply_mass(*(Teuchos::rcp_const_cast<std::vector<Real> >(dual_vec_->getVector())),*vec_);
609  return *dual_vec_;
610  }
611 
612 };
613 
614 template<class Real>
615 class L2VectorDual : public ROL::Vector<Real> {
616 private:
617  Teuchos::RCP<std::vector<Real> > vec_;
618  Teuchos::RCP<BurgersFEM<Real> > fem_;
619 
620  mutable Teuchos::RCP<L2VectorPrimal<Real> > dual_vec_;
621 
622 public:
623  L2VectorDual(const Teuchos::RCP<std::vector<Real> > & vec,
624  const Teuchos::RCP<BurgersFEM<Real> > &fem)
625  : vec_(vec), fem_(fem), dual_vec_(Teuchos::null) {}
626 
627  void set( const ROL::Vector<Real> &x ) {
628  const L2VectorDual &ex = Teuchos::dyn_cast<const L2VectorDual>(x);
629  const std::vector<Real>& xval = *ex.getVector();
630  std::copy(xval.begin(),xval.end(),vec_->begin());
631  }
632 
633  void plus( const ROL::Vector<Real> &x ) {
634  const L2VectorDual &ex = Teuchos::dyn_cast<const L2VectorDual>(x);
635  const std::vector<Real>& xval = *ex.getVector();
636  unsigned dimension = vec_->size();
637  for (unsigned i=0; i<dimension; i++) {
638  (*vec_)[i] += xval[i];
639  }
640  }
641 
642  void scale( const Real alpha ) {
643  unsigned dimension = vec_->size();
644  for (unsigned i=0; i<dimension; i++) {
645  (*vec_)[i] *= alpha;
646  }
647  }
648 
649  Real dot( const ROL::Vector<Real> &x ) const {
650  const L2VectorDual & ex = Teuchos::dyn_cast<const L2VectorDual>(x);
651  const std::vector<Real>& xval = *ex.getVector();
652  unsigned dimension = vec_->size();
653  std::vector<Real> Mx(dimension,0.0);
654  fem_->apply_inverse_mass(Mx,xval);
655  Real val = 0.0;
656  for (unsigned i = 0; i < dimension; i++) {
657  val += Mx[i]*(*vec_)[i];
658  }
659  return val;
660  }
661 
662  Real norm() const {
663  Real val = 0;
664  val = std::sqrt( dot(*this) );
665  return val;
666  }
667 
668  Teuchos::RCP<ROL::Vector<Real> > clone() const {
669  return Teuchos::rcp( new L2VectorDual( Teuchos::rcp(new std::vector<Real>(vec_->size(),0.0)),fem_));
670  }
671 
672  Teuchos::RCP<const std::vector<Real> > getVector() const {
673  return vec_;
674  }
675 
676  Teuchos::RCP<std::vector<Real> > getVector() {
677  return vec_;
678  }
679 
680  Teuchos::RCP<ROL::Vector<Real> > basis( const int i ) const {
681  Teuchos::RCP<L2VectorDual> e
682  = Teuchos::rcp( new L2VectorDual( Teuchos::rcp(new std::vector<Real>(vec_->size(),0.0)),fem_));
683  (*e->getVector())[i] = 1.0;
684  return e;
685  }
686 
687  int dimension() const {
688  return vec_->size();
689  }
690 
691  const ROL::Vector<Real>& dual() const {
692  dual_vec_ = Teuchos::rcp(new L2VectorPrimal<Real>(
693  Teuchos::rcp(new std::vector<Real>(*vec_)),fem_));
694 
695  fem_->apply_inverse_mass(*(Teuchos::rcp_const_cast<std::vector<Real> >(dual_vec_->getVector())),*vec_);
696  return *dual_vec_;
697  }
698 
699 };
700 
701 template<class Real>
702 class H1VectorPrimal : public ROL::Vector<Real> {
703 private:
704  Teuchos::RCP<std::vector<Real> > vec_;
705  Teuchos::RCP<BurgersFEM<Real> > fem_;
706 
707  mutable Teuchos::RCP<H1VectorDual<Real> > dual_vec_;
708 
709 public:
710  H1VectorPrimal(const Teuchos::RCP<std::vector<Real> > & vec,
711  const Teuchos::RCP<BurgersFEM<Real> > &fem)
712  : vec_(vec), fem_(fem), dual_vec_(Teuchos::null) {}
713 
714  void set( const ROL::Vector<Real> &x ) {
715  const H1VectorPrimal &ex = Teuchos::dyn_cast<const H1VectorPrimal>(x);
716  const std::vector<Real>& xval = *ex.getVector();
717  std::copy(xval.begin(),xval.end(),vec_->begin());
718  }
719 
720  void plus( const ROL::Vector<Real> &x ) {
721  const H1VectorPrimal &ex = Teuchos::dyn_cast<const H1VectorPrimal>(x);
722  const std::vector<Real>& xval = *ex.getVector();
723  unsigned dimension = vec_->size();
724  for (unsigned i=0; i<dimension; i++) {
725  (*vec_)[i] += xval[i];
726  }
727  }
728 
729  void scale( const Real alpha ) {
730  unsigned dimension = vec_->size();
731  for (unsigned i=0; i<dimension; i++) {
732  (*vec_)[i] *= alpha;
733  }
734  }
735 
736  Real dot( const ROL::Vector<Real> &x ) const {
737  const H1VectorPrimal & ex = Teuchos::dyn_cast<const H1VectorPrimal>(x);
738  const std::vector<Real>& xval = *ex.getVector();
739  return fem_->compute_H1_dot(xval,*vec_);
740  }
741 
742  Real norm() const {
743  Real val = 0;
744  val = std::sqrt( dot(*this) );
745  return val;
746  }
747 
748  Teuchos::RCP<ROL::Vector<Real> > clone() const {
749  return Teuchos::rcp( new H1VectorPrimal( Teuchos::rcp(new std::vector<Real>(vec_->size(),0.0)),fem_));
750  }
751 
752  Teuchos::RCP<const std::vector<Real> > getVector() const {
753  return vec_;
754  }
755 
756  Teuchos::RCP<std::vector<Real> > getVector() {
757  return vec_;
758  }
759 
760  Teuchos::RCP<ROL::Vector<Real> > basis( const int i ) const {
761  Teuchos::RCP<H1VectorPrimal> e
762  = Teuchos::rcp( new H1VectorPrimal( Teuchos::rcp(new std::vector<Real>(vec_->size(),0.0)),fem_));
763  (*e->getVector())[i] = 1.0;
764  return e;
765  }
766 
767  int dimension() const {
768  return vec_->size();
769  }
770 
771  const ROL::Vector<Real>& dual() const {
772  dual_vec_ = Teuchos::rcp(new H1VectorDual<Real>(
773  Teuchos::rcp(new std::vector<Real>(*vec_)),fem_));
774 
775  fem_->apply_H1(*(Teuchos::rcp_const_cast<std::vector<Real> >(dual_vec_->getVector())),*vec_);
776  return *dual_vec_;
777  }
778 
779 };
780 
781 template<class Real>
782 class H1VectorDual : public ROL::Vector<Real> {
783 private:
784  Teuchos::RCP<std::vector<Real> > vec_;
785  Teuchos::RCP<BurgersFEM<Real> > fem_;
786 
787  mutable Teuchos::RCP<H1VectorPrimal<Real> > dual_vec_;
788 
789 public:
790  H1VectorDual(const Teuchos::RCP<std::vector<Real> > & vec,
791  const Teuchos::RCP<BurgersFEM<Real> > &fem)
792  : vec_(vec), fem_(fem), dual_vec_(Teuchos::null) {}
793 
794  void set( const ROL::Vector<Real> &x ) {
795  const H1VectorDual &ex = Teuchos::dyn_cast<const H1VectorDual>(x);
796  const std::vector<Real>& xval = *ex.getVector();
797  std::copy(xval.begin(),xval.end(),vec_->begin());
798  }
799 
800  void plus( const ROL::Vector<Real> &x ) {
801  const H1VectorDual &ex = Teuchos::dyn_cast<const H1VectorDual>(x);
802  const std::vector<Real>& xval = *ex.getVector();
803  unsigned dimension = vec_->size();
804  for (unsigned i=0; i<dimension; i++) {
805  (*vec_)[i] += xval[i];
806  }
807  }
808 
809  void scale( const Real alpha ) {
810  unsigned dimension = vec_->size();
811  for (unsigned i=0; i<dimension; i++) {
812  (*vec_)[i] *= alpha;
813  }
814  }
815 
816  Real dot( const ROL::Vector<Real> &x ) const {
817  const H1VectorDual & ex = Teuchos::dyn_cast<const H1VectorDual>(x);
818  const std::vector<Real>& xval = *ex.getVector();
819  unsigned dimension = vec_->size();
820  std::vector<Real> Mx(dimension,0.0);
821  fem_->apply_inverse_H1(Mx,xval);
822  Real val = 0.0;
823  for (unsigned i = 0; i < dimension; i++) {
824  val += Mx[i]*(*vec_)[i];
825  }
826  return val;
827  }
828 
829  Real norm() const {
830  Real val = 0;
831  val = std::sqrt( dot(*this) );
832  return val;
833  }
834 
835  Teuchos::RCP<ROL::Vector<Real> > clone() const {
836  return Teuchos::rcp( new H1VectorDual( Teuchos::rcp(new std::vector<Real>(vec_->size(),0.0)),fem_));
837  }
838 
839  Teuchos::RCP<const std::vector<Real> > getVector() const {
840  return vec_;
841  }
842 
843  Teuchos::RCP<std::vector<Real> > getVector() {
844  return vec_;
845  }
846 
847  Teuchos::RCP<ROL::Vector<Real> > basis( const int i ) const {
848  Teuchos::RCP<H1VectorDual> e
849  = Teuchos::rcp( new H1VectorDual( Teuchos::rcp(new std::vector<Real>(vec_->size(),0.0)),fem_));
850  (*e->getVector())[i] = 1.0;
851  return e;
852  }
853 
854  int dimension() const {
855  return vec_->size();
856  }
857 
858  const ROL::Vector<Real>& dual() const {
859  dual_vec_ = Teuchos::rcp(new H1VectorPrimal<Real>(
860  Teuchos::rcp(new std::vector<Real>(*vec_)),fem_));
861 
862  fem_->apply_inverse_H1(*(Teuchos::rcp_const_cast<std::vector<Real> >(dual_vec_->getVector())),*vec_);
863  return *dual_vec_;
864  }
865 
866 };
867 
868 template<class Real>
870 private:
871 
874 
877 
880 
881  typedef typename std::vector<Real>::size_type uint;
882 
883  Teuchos::RCP<BurgersFEM<Real> > fem_;
884  bool useHessian_;
885 
886 public:
887  EqualityConstraint_BurgersControl(Teuchos::RCP<BurgersFEM<Real> > &fem, bool useHessian = true)
888  : fem_(fem), useHessian_(useHessian) {}
889 
891  const ROL::Vector<Real> &z, Real &tol) {
892  Teuchos::RCP<std::vector<Real> > cp =
893  Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<PrimalConstraintVector>(c)).getVector());
894  Teuchos::RCP<const std::vector<Real> > up =
895  (Teuchos::dyn_cast<PrimalStateVector>(const_cast<ROL::Vector<Real> &>(u))).getVector();
896  Teuchos::RCP<const std::vector<Real> > zp =
897  (Teuchos::dyn_cast<PrimalControlVector>(const_cast<ROL::Vector<Real> &>(z))).getVector();
898 
899  const std::vector<Real> param
901  fem_->set_problem_data(param[0],param[1],param[2],param[3]);
902 
903  fem_->compute_residual(*cp,*up,*zp);
904  }
905 
907  const ROL::Vector<Real> &z, Real &tol) {
908  Teuchos::RCP<std::vector<Real> > jvp =
909  Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<PrimalConstraintVector>(jv)).getVector());
910  Teuchos::RCP<const std::vector<Real> > vp =
911  (Teuchos::dyn_cast<PrimalStateVector>(const_cast<ROL::Vector<Real> &>(v))).getVector();
912  Teuchos::RCP<const std::vector<Real> > up =
913  (Teuchos::dyn_cast<PrimalStateVector>(const_cast<ROL::Vector<Real> &>(u))).getVector();
914  Teuchos::RCP<const std::vector<Real> > zp =
915  (Teuchos::dyn_cast<PrimalControlVector>(const_cast<ROL::Vector<Real> &>(z))).getVector();
916 
917  const std::vector<Real> param
919  fem_->set_problem_data(param[0],param[1],param[2],param[3]);
920 
921  fem_->apply_pde_jacobian(*jvp,*vp,*up,*zp);
922  }
923 
925  const ROL::Vector<Real> &z, Real &tol) {
926  Teuchos::RCP<std::vector<Real> > jvp =
927  Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<PrimalConstraintVector>(jv)).getVector());
928  Teuchos::RCP<const std::vector<Real> > vp =
929  (Teuchos::dyn_cast<PrimalControlVector>(const_cast<ROL::Vector<Real> &>(v))).getVector();
930  Teuchos::RCP<const std::vector<Real> > up =
931  (Teuchos::dyn_cast<PrimalStateVector>(const_cast<ROL::Vector<Real> &>(u))).getVector();
932  Teuchos::RCP<const std::vector<Real> > zp =
933  (Teuchos::dyn_cast<PrimalControlVector>(const_cast<ROL::Vector<Real> &>(z))).getVector();
934 
935  const std::vector<Real> param
937  fem_->set_problem_data(param[0],param[1],param[2],param[3]);
938 
939  fem_->apply_control_jacobian(*jvp,*vp,*up,*zp);
940  }
941 
943  const ROL::Vector<Real> &z, Real &tol) {
944  Teuchos::RCP<std::vector<Real> > ijvp =
945  Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<PrimalStateVector>(ijv)).getVector());
946  Teuchos::RCP<const std::vector<Real> > vp =
947  (Teuchos::dyn_cast<PrimalConstraintVector>(const_cast<ROL::Vector<Real> &>(v))).getVector();
948  Teuchos::RCP<const std::vector<Real> > up =
949  (Teuchos::dyn_cast<PrimalStateVector>(const_cast<ROL::Vector<Real> &>(u))).getVector();
950  Teuchos::RCP<const std::vector<Real> > zp =
951  (Teuchos::dyn_cast<PrimalControlVector>(const_cast<ROL::Vector<Real> &>(z))).getVector();
952 
953  const std::vector<Real> param
955  fem_->set_problem_data(param[0],param[1],param[2],param[3]);
956 
957  fem_->apply_inverse_pde_jacobian(*ijvp,*vp,*up,*zp);
958  }
959 
961  const ROL::Vector<Real> &z, Real &tol) {
962  Teuchos::RCP<std::vector<Real> > jvp =
963  Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<DualStateVector>(ajv)).getVector());
964  Teuchos::RCP<const std::vector<Real> > vp =
965  (Teuchos::dyn_cast<DualConstraintVector>(const_cast<ROL::Vector<Real> &>(v))).getVector();
966  Teuchos::RCP<const std::vector<Real> > up =
967  (Teuchos::dyn_cast<PrimalStateVector>(const_cast<ROL::Vector<Real> &>(u))).getVector();
968  Teuchos::RCP<const std::vector<Real> > zp =
969  (Teuchos::dyn_cast<PrimalControlVector>(const_cast<ROL::Vector<Real> &>(z))).getVector();
970 
971  const std::vector<Real> param
973  fem_->set_problem_data(param[0],param[1],param[2],param[3]);
974 
975  fem_->apply_adjoint_pde_jacobian(*jvp,*vp,*up,*zp);
976  }
977 
979  const ROL::Vector<Real> &z, Real &tol) {
980  Teuchos::RCP<std::vector<Real> > jvp =
981  Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<DualControlVector>(jv)).getVector());
982  Teuchos::RCP<const std::vector<Real> > vp =
983  (Teuchos::dyn_cast<DualConstraintVector>(const_cast<ROL::Vector<Real> &>(v))).getVector();
984  Teuchos::RCP<const std::vector<Real> > up =
985  (Teuchos::dyn_cast<PrimalStateVector>(const_cast<ROL::Vector<Real> &>(u))).getVector();
986  Teuchos::RCP<const std::vector<Real> > zp =
987  (Teuchos::dyn_cast<PrimalControlVector>(const_cast<ROL::Vector<Real> &>(z))).getVector();
988 
989  const std::vector<Real> param
991  fem_->set_problem_data(param[0],param[1],param[2],param[3]);
992 
993  fem_->apply_adjoint_control_jacobian(*jvp,*vp,*up,*zp);
994  }
995 
997  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
998  Teuchos::RCP<std::vector<Real> > iajvp =
999  Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<DualConstraintVector>(iajv)).getVector());
1000  Teuchos::RCP<const std::vector<Real> > vp =
1001  (Teuchos::dyn_cast<DualStateVector>(const_cast<ROL::Vector<Real> &>(v))).getVector();
1002  Teuchos::RCP<const std::vector<Real> > up =
1003  (Teuchos::dyn_cast<PrimalStateVector>(const_cast<ROL::Vector<Real> &>(u))).getVector();
1004  Teuchos::RCP<const std::vector<Real> > zp =
1005  (Teuchos::dyn_cast<PrimalControlVector>(const_cast<ROL::Vector<Real> &>(z))).getVector();
1006 
1007  const std::vector<Real> param
1009  fem_->set_problem_data(param[0],param[1],param[2],param[3]);
1010 
1011  fem_->apply_inverse_adjoint_pde_jacobian(*iajvp,*vp,*up,*zp);
1012  }
1013 
1015  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
1016  if ( useHessian_ ) {
1017  Teuchos::RCP<std::vector<Real> > ahwvp =
1018  Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<DualStateVector>(ahwv)).getVector());
1019  Teuchos::RCP<const std::vector<Real> > wp =
1020  (Teuchos::dyn_cast<DualConstraintVector>(const_cast<ROL::Vector<Real> &>(w))).getVector();
1021  Teuchos::RCP<const std::vector<Real> > vp =
1022  (Teuchos::dyn_cast<PrimalStateVector>(const_cast<ROL::Vector<Real> &>(v))).getVector();
1023  Teuchos::RCP<const std::vector<Real> > up =
1024  (Teuchos::dyn_cast<PrimalStateVector>(const_cast<ROL::Vector<Real> &>(u))).getVector();
1025  Teuchos::RCP<const std::vector<Real> > zp =
1026  (Teuchos::dyn_cast<PrimalControlVector>(const_cast<ROL::Vector<Real> &>(z))).getVector();
1027 
1028  const std::vector<Real> param
1030  fem_->set_problem_data(param[0],param[1],param[2],param[3]);
1031 
1032  fem_->apply_adjoint_pde_hessian(*ahwvp,*wp,*vp,*up,*zp);
1033  }
1034  else {
1035  ahwv.zero();
1036  }
1037  }
1038 
1040  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
1041  if ( useHessian_ ) {
1042  Teuchos::RCP<std::vector<Real> > ahwvp =
1043  Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<DualControlVector>(ahwv)).getVector());
1044  Teuchos::RCP<const std::vector<Real> > wp =
1045  (Teuchos::dyn_cast<DualConstraintVector>(const_cast<ROL::Vector<Real> &>(w))).getVector();
1046  Teuchos::RCP<const std::vector<Real> > vp =
1047  (Teuchos::dyn_cast<PrimalStateVector>(const_cast<ROL::Vector<Real> &>(v))).getVector();
1048  Teuchos::RCP<const std::vector<Real> > up =
1049  (Teuchos::dyn_cast<PrimalStateVector>(const_cast<ROL::Vector<Real> &>(u))).getVector();
1050  Teuchos::RCP<const std::vector<Real> > zp =
1051  (Teuchos::dyn_cast<PrimalControlVector>(const_cast<ROL::Vector<Real> &>(z))).getVector();
1052 
1053  const std::vector<Real> param
1055  fem_->set_problem_data(param[0],param[1],param[2],param[3]);
1056 
1057  fem_->apply_adjoint_control_pde_hessian(*ahwvp,*wp,*vp,*up,*zp);
1058  }
1059  else {
1060  ahwv.zero();
1061  }
1062  }
1064  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
1065  if ( useHessian_ ) {
1066  Teuchos::RCP<std::vector<Real> > ahwvp =
1067  Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<DualStateVector>(ahwv)).getVector());
1068  Teuchos::RCP<const std::vector<Real> > wp =
1069  (Teuchos::dyn_cast<DualConstraintVector>(const_cast<ROL::Vector<Real> &>(w))).getVector();
1070  Teuchos::RCP<const std::vector<Real> > vp =
1071  (Teuchos::dyn_cast<PrimalControlVector>(const_cast<ROL::Vector<Real> &>(v))).getVector();
1072  Teuchos::RCP<const std::vector<Real> > up =
1073  (Teuchos::dyn_cast<PrimalStateVector>(const_cast<ROL::Vector<Real> &>(u))).getVector();
1074  Teuchos::RCP<const std::vector<Real> > zp =
1075  (Teuchos::dyn_cast<PrimalControlVector>(const_cast<ROL::Vector<Real> &>(z))).getVector();
1076 
1077  const std::vector<Real> param
1079  fem_->set_problem_data(param[0],param[1],param[2],param[3]);
1080 
1081  fem_->apply_adjoint_pde_control_hessian(*ahwvp,*wp,*vp,*up,*zp);
1082  }
1083  else {
1084  ahwv.zero();
1085  }
1086  }
1088  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
1089  if ( useHessian_ ) {
1090  Teuchos::RCP<std::vector<Real> > ahwvp =
1091  Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<DualControlVector>(ahwv)).getVector());
1092  Teuchos::RCP<const std::vector<Real> > wp =
1093  (Teuchos::dyn_cast<DualConstraintVector>(const_cast<ROL::Vector<Real> &>(w))).getVector();
1094  Teuchos::RCP<const std::vector<Real> > vp =
1095  (Teuchos::dyn_cast<PrimalControlVector>(const_cast<ROL::Vector<Real> &>(v))).getVector();
1096  Teuchos::RCP<const std::vector<Real> > up =
1097  (Teuchos::dyn_cast<PrimalStateVector>(const_cast<ROL::Vector<Real> &>(u))).getVector();
1098  Teuchos::RCP<const std::vector<Real> > zp =
1099  (Teuchos::dyn_cast<PrimalControlVector>(const_cast<ROL::Vector<Real> &>(z))).getVector();
1100 
1101  const std::vector<Real> param
1103  fem_->set_problem_data(param[0],param[1],param[2],param[3]);
1104 
1105  fem_->apply_adjoint_control_hessian(*ahwvp,*wp,*vp,*up,*zp);
1106  }
1107  else {
1108  ahwv.zero();
1109  }
1110  }
1111 };
1112 
1113 template<class Real>
1114 class Objective_BurgersControl : public ROL::Objective_SimOpt<Real> {
1115 private:
1116 
1119 
1122 
1123  Teuchos::RCP<BurgersFEM<Real> > fem_;
1124 
1125  Real x_;
1126  std::vector<int> indices_;
1127 
1128 public:
1129  Objective_BurgersControl(const Teuchos::RCP<BurgersFEM<Real> > &fem,
1130  Real x = 0.0) : fem_(fem), x_(x) {
1131  for (int i = 1; i < fem_->num_dof()+1; i++) {
1132  if ( (Real)i*(fem_->mesh_spacing()) >= x_ ) {
1133  indices_.push_back(i-1);
1134  }
1135  }
1136  }
1137 
1139 
1140  Real value( const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
1141  Teuchos::RCP<const std::vector<Real> > up =
1142  (Teuchos::dyn_cast<PrimalStateVector>(const_cast<ROL::Vector<Real> &>(u))).getVector();
1143 
1144 // const std::vector<Real> param
1145 // = ROL::Objective_SimOpt<Real>::getParameter();
1146 // fem_->set_problem_data(param[0],param[1],param[2],param[3]);
1147 // Real nu = fem_->get_viscosity();
1148 //
1149 // return 0.5*nu*fem_->compute_H1_dot(*up,*up);
1150 
1151  Real val = 0.5*((((Real)indices_[0]+1.)*(fem_->mesh_spacing())-x_)
1152  *(x_+(2.-((Real)indices_[0]+1.))*(fem_->mesh_spacing()))/(fem_->mesh_spacing())
1153  +(fem_->mesh_spacing())) * (*up)[indices_[0]];
1154  for (uint i = 1; i < indices_.size(); i++) {
1155  val += (fem_->mesh_spacing())*(*up)[indices_[i]];
1156  }
1157  return -val;
1158  }
1159 
1160  void gradient_1( ROL::Vector<Real> &g, const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
1161  Teuchos::RCP<std::vector<Real> > gp =
1162  Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<DualStateVector>(g)).getVector());
1163  Teuchos::RCP<const std::vector<Real> > up =
1164  (Teuchos::dyn_cast<PrimalStateVector>(const_cast<ROL::Vector<Real> &>(u))).getVector();
1165 
1166 // const std::vector<Real> param
1167 // = ROL::Objective_SimOpt<Real>::getParameter();
1168 // fem_->set_problem_data(param[0],param[1],param[2],param[3]);
1169 // Real nu = fem_->get_viscosity();
1170 //
1171 // fem_->apply_H1(*gp,*up);
1172 // g.scale(nu);
1173 
1174  g.zero();
1175  (*gp)[indices_[0]] = -0.5*((((Real)indices_[0]+1.)*(fem_->mesh_spacing())-x_)
1176  *(x_+(2.-((Real)indices_[0]+1.))*(fem_->mesh_spacing()))/(fem_->mesh_spacing())
1177  +(fem_->mesh_spacing()));
1178 
1179 
1180  for (uint i = 1; i < indices_.size(); i++) {
1181  (*gp)[indices_[i]] = -(fem_->mesh_spacing());
1182  }
1183  }
1184 
1185  void gradient_2( ROL::Vector<Real> &g, const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
1186  g.zero();
1187  }
1188 
1190  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
1191 // Teuchos::RCP<std::vector<Real> > hvp =
1192 // Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<DualStateVector>(hv)).getVector());
1193 // Teuchos::RCP<const std::vector<Real> > vp =
1194 // (Teuchos::dyn_cast<PrimalStateVector>(const_cast<ROL::Vector<Real> &>(v))).getVector();
1195 //
1196 // const std::vector<Real> param
1197 // = ROL::Objective_SimOpt<Real>::getParameter();
1198 // fem_->set_problem_data(param[0],param[1],param[2],param[3]);
1199 // Real nu = fem_->get_viscosity();
1200 //
1201 // fem_->apply_H1(*hvp,*vp);
1202 // hv.scale(nu);
1203 
1204  hv.zero();
1205  }
1206 
1208  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
1209  hv.zero();
1210  }
1211 
1213  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
1214  hv.zero();
1215  }
1216 
1218  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
1219  hv.zero();
1220  }
1221 };
1222 
1223 template<class Real>
1224 class L2BoundConstraint : public ROL::BoundConstraint<Real> {
1225 private:
1226  int dim_;
1227  std::vector<Real> x_lo_;
1228  std::vector<Real> x_up_;
1229  Real min_diff_;
1230  Real scale_;
1231  Teuchos::RCP<BurgersFEM<Real> > fem_;
1232 
1233  void cast_vector(Teuchos::RCP<std::vector<Real> > &xvec,
1234  ROL::Vector<Real> &x) const {
1235  try {
1236  xvec = Teuchos::rcp_const_cast<std::vector<Real> >(
1237  (Teuchos::dyn_cast<L2VectorPrimal<Real> >(x)).getVector());
1238  }
1239  catch (std::exception &e) {
1240  xvec = Teuchos::rcp_const_cast<std::vector<Real> >(
1241  (Teuchos::dyn_cast<L2VectorDual<Real> >(x)).getVector());
1242  }
1243  }
1244 
1245  void cast_const_vector(Teuchos::RCP<const std::vector<Real> > &xvec,
1246  const ROL::Vector<Real> &x) const {
1247  try {
1248  xvec = (Teuchos::dyn_cast<L2VectorPrimal<Real> >(
1249  const_cast<ROL::Vector<Real> &>(x))).getVector();
1250  }
1251  catch (std::exception &e) {
1252  xvec = (Teuchos::dyn_cast<L2VectorDual<Real> >(
1253  const_cast<ROL::Vector<Real> &>(x))).getVector();
1254  }
1255  }
1256 
1257  void axpy(std::vector<Real> &out, const Real a,
1258  const std::vector<Real> &x, const std::vector<Real> &y) const{
1259  out.resize(dim_,0.0);
1260  for (unsigned i = 0; i < dim_; i++) {
1261  out[i] = a*x[i] + y[i];
1262  }
1263  }
1264 
1265  void projection(std::vector<Real> &x) {
1266  for ( int i = 0; i < dim_; i++ ) {
1267  x[i] = std::max(x_lo_[i],std::min(x_up_[i],x[i]));
1268  }
1269  }
1270 
1271 public:
1272  L2BoundConstraint(std::vector<Real> &l, std::vector<Real> &u,
1273  const Teuchos::RCP<BurgersFEM<Real> > &fem, Real scale = 1.0)
1274  : x_lo_(l), x_up_(u), scale_(scale), fem_(fem) {
1275  dim_ = x_lo_.size();
1276  for ( int i = 0; i < dim_; i++ ) {
1277  if ( i == 0 ) {
1278  min_diff_ = x_up_[i] - x_lo_[i];
1279  }
1280  else {
1281  min_diff_ = ( (min_diff_ < (x_up_[i] - x_lo_[i])) ? min_diff_ : (x_up_[i] - x_lo_[i]) );
1282  }
1283  }
1284  min_diff_ *= 0.5;
1285  }
1286 
1287  bool isFeasible( const ROL::Vector<Real> &x ) {
1288  Teuchos::RCP<const std::vector<Real> > ex; cast_const_vector(ex,x);
1289  bool val = true;
1290  int cnt = 1;
1291  for ( int i = 0; i < dim_; i++ ) {
1292  if ( (*ex)[i] >= x_lo_[i] && (*ex)[i] <= x_up_[i] ) { cnt *= 1; }
1293  else { cnt *= 0; }
1294  }
1295  if ( cnt == 0 ) { val = false; }
1296  return val;
1297  }
1298 
1300  Teuchos::RCP<std::vector<Real> > ex; cast_vector(ex,x);
1301  projection(*ex);
1302  }
1303 
1305  Teuchos::RCP<const std::vector<Real> > ex; cast_const_vector(ex,x);
1306  Teuchos::RCP<std::vector<Real> > ev; cast_vector(ev,v);
1307  Real epsn = std::min(scale_*eps,min_diff_);
1308  for ( int i = 0; i < dim_; i++ ) {
1309  if ( ((*ex)[i] <= x_lo_[i]+epsn) ) {
1310  (*ev)[i] = 0.0;
1311  }
1312  }
1313  }
1314 
1316  Teuchos::RCP<const std::vector<Real> > ex; cast_const_vector(ex,x);
1317  Teuchos::RCP<std::vector<Real> > ev; cast_vector(ev,v);
1318  Real epsn = std::min(scale_*eps,min_diff_);
1319  for ( int i = 0; i < dim_; i++ ) {
1320  if ( ((*ex)[i] >= x_up_[i]-epsn) ) {
1321  (*ev)[i] = 0.0;
1322  }
1323  }
1324  }
1325 
1326  void pruneActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &x, Real eps) {
1327  Teuchos::RCP<const std::vector<Real> > ex; cast_const_vector(ex,x);
1328  Teuchos::RCP<std::vector<Real> > ev; cast_vector(ev,v);
1329  Real epsn = std::min(scale_*eps,min_diff_);
1330  for ( int i = 0; i < dim_; i++ ) {
1331  if ( ((*ex)[i] <= x_lo_[i]+epsn) ||
1332  ((*ex)[i] >= x_up_[i]-epsn) ) {
1333  (*ev)[i] = 0.0;
1334  }
1335  }
1336  }
1337 
1339  Teuchos::RCP<const std::vector<Real> > ex; cast_const_vector(ex,x);
1340  Teuchos::RCP<const std::vector<Real> > eg; cast_const_vector(eg,g);
1341  Teuchos::RCP<std::vector<Real> > ev; cast_vector(ev,v);
1342  Real epsn = std::min(scale_*eps,min_diff_);
1343  for ( int i = 0; i < dim_; i++ ) {
1344  if ( ((*ex)[i] <= x_lo_[i]+epsn && (*eg)[i] > 0.0) ) {
1345  (*ev)[i] = 0.0;
1346  }
1347  }
1348  }
1349 
1351  Teuchos::RCP<const std::vector<Real> > ex; cast_const_vector(ex,x);
1352  Teuchos::RCP<const std::vector<Real> > eg; cast_const_vector(eg,g);
1353  Teuchos::RCP<std::vector<Real> > ev; cast_vector(ev,v);
1354  Real epsn = std::min(scale_*eps,min_diff_);
1355  for ( int i = 0; i < dim_; i++ ) {
1356  if ( ((*ex)[i] >= x_up_[i]-epsn && (*eg)[i] < 0.0) ) {
1357  (*ev)[i] = 0.0;
1358  }
1359  }
1360  }
1361 
1362  void pruneActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &g, const ROL::Vector<Real> &x, Real eps) {
1363  Teuchos::RCP<const std::vector<Real> > ex; cast_const_vector(ex,x);
1364  Teuchos::RCP<const std::vector<Real> > eg; cast_const_vector(eg,g);
1365  Teuchos::RCP<std::vector<Real> > ev; cast_vector(ev,v);
1366  Real epsn = std::min(scale_*eps,min_diff_);
1367  for ( int i = 0; i < dim_; i++ ) {
1368  if ( ((*ex)[i] <= x_lo_[i]+epsn && (*eg)[i] > 0.0) ||
1369  ((*ex)[i] >= x_up_[i]-epsn && (*eg)[i] < 0.0) ) {
1370  (*ev)[i] = 0.0;
1371  }
1372  }
1373  }
1374 
1376  Teuchos::RCP<std::vector<Real> > us = Teuchos::rcp( new std::vector<Real>(dim_,0.0) );
1377  us->assign(x_up_.begin(),x_up_.end());
1378  Teuchos::RCP<ROL::Vector<Real> > up = Teuchos::rcp( new L2VectorPrimal<Real>(us,fem_) );
1379  u.set(*up);
1380  }
1381 
1383  Teuchos::RCP<std::vector<Real> > ls = Teuchos::rcp( new std::vector<Real>(dim_,0.0) );
1384  ls->assign(x_lo_.begin(),x_lo_.end());
1385  Teuchos::RCP<ROL::Vector<Real> > lp = Teuchos::rcp( new L2VectorPrimal<Real>(ls,fem_) );
1386  l.set(*lp);
1387  }
1388 };
1389 
1390 template<class Real>
1391 class H1BoundConstraint : public ROL::BoundConstraint<Real> {
1392 private:
1393  int dim_;
1394  std::vector<Real> x_lo_;
1395  std::vector<Real> x_up_;
1396  Real min_diff_;
1397  Real scale_;
1398  Teuchos::RCP<BurgersFEM<Real> > fem_;
1399 
1400  void cast_vector(Teuchos::RCP<std::vector<Real> > &xvec,
1401  ROL::Vector<Real> &x) const {
1402  try {
1403  xvec = Teuchos::rcp_const_cast<std::vector<Real> >(
1404  (Teuchos::dyn_cast<H1VectorPrimal<Real> >(x)).getVector());
1405  }
1406  catch (std::exception &e) {
1407  xvec = Teuchos::rcp_const_cast<std::vector<Real> >(
1408  (Teuchos::dyn_cast<H1VectorDual<Real> >(x)).getVector());
1409  }
1410  }
1411 
1412  void cast_const_vector(Teuchos::RCP<const std::vector<Real> > &xvec,
1413  const ROL::Vector<Real> &x) const {
1414  try {
1415  xvec = (Teuchos::dyn_cast<H1VectorPrimal<Real> >(
1416  const_cast<ROL::Vector<Real> &>(x))).getVector();
1417  }
1418  catch (std::exception &e) {
1419  xvec = (Teuchos::dyn_cast<H1VectorDual<Real> >(
1420  const_cast<ROL::Vector<Real> &>(x))).getVector();
1421  }
1422  }
1423 
1424  void axpy(std::vector<Real> &out, const Real a,
1425  const std::vector<Real> &x, const std::vector<Real> &y) const{
1426  out.resize(dim_,0.0);
1427  for (unsigned i = 0; i < dim_; i++) {
1428  out[i] = a*x[i] + y[i];
1429  }
1430  }
1431 
1432  void projection(std::vector<Real> &x) {
1433  for ( int i = 0; i < dim_; i++ ) {
1434  x[i] = std::max(x_lo_[i],std::min(x_up_[i],x[i]));
1435  }
1436  }
1437 
1438 public:
1439  H1BoundConstraint(std::vector<Real> &l, std::vector<Real> &u,
1440  const Teuchos::RCP<BurgersFEM<Real> > &fem, Real scale = 1.0)
1441  : x_lo_(l), x_up_(u), scale_(scale), fem_(fem) {
1442  dim_ = x_lo_.size();
1443  for ( int i = 0; i < dim_; i++ ) {
1444  if ( i == 0 ) {
1445  min_diff_ = x_up_[i] - x_lo_[i];
1446  }
1447  else {
1448  min_diff_ = ( (min_diff_ < (x_up_[i] - x_lo_[i])) ? min_diff_ : (x_up_[i] - x_lo_[i]) );
1449  }
1450  }
1451  min_diff_ *= 0.5;
1452  }
1453 
1454  bool isFeasible( const ROL::Vector<Real> &x ) {
1455  Teuchos::RCP<const std::vector<Real> > ex; cast_const_vector(ex,x);
1456  bool val = true;
1457  int cnt = 1;
1458  for ( int i = 0; i < dim_; i++ ) {
1459  if ( (*ex)[i] >= x_lo_[i] && (*ex)[i] <= x_up_[i] ) { cnt *= 1; }
1460  else { cnt *= 0; }
1461  }
1462  if ( cnt == 0 ) { val = false; }
1463  return val;
1464  }
1465 
1467  Teuchos::RCP<std::vector<Real> > ex; cast_vector(ex,x);
1468  projection(*ex);
1469  }
1470 
1472  Teuchos::RCP<const std::vector<Real> > ex; cast_const_vector(ex,x);
1473  Teuchos::RCP<std::vector<Real> > ev; cast_vector(ev,v);
1474  Real epsn = std::min(scale_*eps,min_diff_);
1475  for ( int i = 0; i < dim_; i++ ) {
1476  if ( ((*ex)[i] <= x_lo_[i]+epsn) ) {
1477  (*ev)[i] = 0.0;
1478  }
1479  }
1480  }
1481 
1483  Teuchos::RCP<const std::vector<Real> > ex; cast_const_vector(ex,x);
1484  Teuchos::RCP<std::vector<Real> > ev; cast_vector(ev,v);
1485  Real epsn = std::min(scale_*eps,min_diff_);
1486  for ( int i = 0; i < dim_; i++ ) {
1487  if ( ((*ex)[i] >= x_up_[i]-epsn) ) {
1488  (*ev)[i] = 0.0;
1489  }
1490  }
1491  }
1492 
1493  void pruneActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &x, Real eps) {
1494  Teuchos::RCP<const std::vector<Real> > ex; cast_const_vector(ex,x);
1495  Teuchos::RCP<std::vector<Real> > ev; cast_vector(ev,v);
1496  Real epsn = std::min(scale_*eps,min_diff_);
1497  for ( int i = 0; i < dim_; i++ ) {
1498  if ( ((*ex)[i] <= x_lo_[i]+epsn) ||
1499  ((*ex)[i] >= x_up_[i]-epsn) ) {
1500  (*ev)[i] = 0.0;
1501  }
1502  }
1503  }
1504 
1506  Teuchos::RCP<const std::vector<Real> > ex; cast_const_vector(ex,x);
1507  Teuchos::RCP<const std::vector<Real> > eg; cast_const_vector(eg,g);
1508  Teuchos::RCP<std::vector<Real> > ev; cast_vector(ev,v);
1509  Real epsn = std::min(scale_*eps,min_diff_);
1510  for ( int i = 0; i < dim_; i++ ) {
1511  if ( ((*ex)[i] <= x_lo_[i]+epsn && (*eg)[i] > 0.0) ) {
1512  (*ev)[i] = 0.0;
1513  }
1514  }
1515  }
1516 
1518  Teuchos::RCP<const std::vector<Real> > ex; cast_const_vector(ex,x);
1519  Teuchos::RCP<const std::vector<Real> > eg; cast_const_vector(eg,g);
1520  Teuchos::RCP<std::vector<Real> > ev; cast_vector(ev,v);
1521  Real epsn = std::min(scale_*eps,min_diff_);
1522  for ( int i = 0; i < dim_; i++ ) {
1523  if ( ((*ex)[i] >= x_up_[i]-epsn && (*eg)[i] < 0.0) ) {
1524  (*ev)[i] = 0.0;
1525  }
1526  }
1527  }
1528 
1529  void pruneActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &g, const ROL::Vector<Real> &x, Real eps) {
1530  Teuchos::RCP<const std::vector<Real> > ex; cast_const_vector(ex,x);
1531  Teuchos::RCP<const std::vector<Real> > eg; cast_const_vector(eg,g);
1532  Teuchos::RCP<std::vector<Real> > ev; cast_vector(ev,v);
1533  Real epsn = std::min(scale_*eps,min_diff_);
1534  for ( int i = 0; i < dim_; i++ ) {
1535  if ( ((*ex)[i] <= x_lo_[i]+epsn && (*eg)[i] > 0.0) ||
1536  ((*ex)[i] >= x_up_[i]-epsn && (*eg)[i] < 0.0) ) {
1537  (*ev)[i] = 0.0;
1538  }
1539  }
1540  }
1541 
1543  Teuchos::RCP<std::vector<Real> > us = Teuchos::rcp( new std::vector<Real>(dim_,0.0) );
1544  us->assign(x_up_.begin(),x_up_.end());
1545  Teuchos::RCP<ROL::Vector<Real> > up = Teuchos::rcp( new H1VectorPrimal<Real>(us,fem_) );
1546  u.set(*up);
1547  }
1548 
1550  Teuchos::RCP<std::vector<Real> > ls = Teuchos::rcp( new std::vector<Real>(dim_,0.0) );
1551  ls->assign(x_lo_.begin(),x_lo_.end());
1552  Teuchos::RCP<ROL::Vector<Real> > lp = Teuchos::rcp( new H1VectorPrimal<Real>(ls,fem_) );
1553  l.set(*lp);
1554  }
1555 };
1556 
1557 template<class Real, class Ordinal>
1558 class L2VectorBatchManager : public ROL::TeuchosBatchManager<Real,Ordinal> {
1559 private:
1560  void cast_vector(Teuchos::RCP<std::vector<Real> > &xvec,
1561  ROL::Vector<Real> &x) const {
1562  try {
1563  xvec = Teuchos::rcp_const_cast<std::vector<Real> >(
1564  (Teuchos::dyn_cast<L2VectorPrimal<Real> >(x)).getVector());
1565  }
1566  catch (std::exception &e) {
1567  xvec = Teuchos::rcp_const_cast<std::vector<Real> >(
1568  (Teuchos::dyn_cast<L2VectorDual<Real> >(x)).getVector());
1569  }
1570  }
1571 
1572 public:
1573  L2VectorBatchManager(const Teuchos::RCP<const Teuchos::Comm<Ordinal> > &comm)
1574  : ROL::TeuchosBatchManager<Real,Ordinal>(comm) {}
1576  Teuchos::RCP<std::vector<Real> > input_ptr;
1577  cast_vector(input_ptr,input);
1578  int dim_i = input_ptr->size();
1579  Teuchos::RCP<std::vector<Real> > output_ptr;
1580  cast_vector(output_ptr,output);
1581  int dim_o = output_ptr->size();
1582  if ( dim_i != dim_o ) {
1583  std::cout << "L2VectorBatchManager: DIMENSION MISMATCH ON RANK "
1585  }
1586  else {
1587  ROL::TeuchosBatchManager<Real,Ordinal>::sumAll(&(*input_ptr)[0],&(*output_ptr)[0],dim_i);
1588  }
1589  }
1590 };
1591 
1592 template<class Real, class Ordinal>
1593 class H1VectorBatchManager : public ROL::TeuchosBatchManager<Real,Ordinal> {
1594 private:
1595  void cast_vector(Teuchos::RCP<std::vector<Real> > &xvec,
1596  ROL::Vector<Real> &x) const {
1597  try {
1598  xvec = Teuchos::rcp_const_cast<std::vector<Real> >(
1599  (Teuchos::dyn_cast<H1VectorPrimal<Real> >(x)).getVector());
1600  }
1601  catch (std::exception &e) {
1602  xvec = Teuchos::rcp_const_cast<std::vector<Real> >(
1603  (Teuchos::dyn_cast<H1VectorDual<Real> >(x)).getVector());
1604  }
1605  }
1606 
1607 public:
1608  H1VectorBatchManager(const Teuchos::RCP<const Teuchos::Comm<Ordinal> > &comm)
1609  : ROL::TeuchosBatchManager<Real,Ordinal>(comm) {}
1611  Teuchos::RCP<std::vector<Real> > input_ptr;
1612  cast_vector(input_ptr,input);
1613  int dim_i = input_ptr->size();
1614  Teuchos::RCP<std::vector<Real> > output_ptr;
1615  cast_vector(output_ptr,output);
1616  int dim_o = output_ptr->size();
1617  if ( dim_i != dim_o ) {
1618  std::cout << "H1VectorBatchManager: DIMENSION MISMATCH ON RANK "
1620  }
1621  else {
1622  ROL::TeuchosBatchManager<Real,Ordinal>::sumAll(&(*input_ptr)[0],&(*output_ptr)[0],dim_i);
1623  }
1624  }
1625 };
1626 
1627 template<class Real>
1628 Real random(const Teuchos::RCP<const Teuchos::Comm<int> > &comm) {
1629  Real val = 0.0;
1630  if ( Teuchos::rank<int>(*comm)==0 ) {
1631  val = (Real)rand()/(Real)RAND_MAX;
1632  }
1633  Teuchos::broadcast<int,Real>(*comm,0,1,&val);
1634  return val;
1635 }
BurgersFEM(int nx=128, Real nl=1.0, Real cH1=1.0, Real cL2=1.0)
Definition: example_07.hpp:133
Provides the interface to evaluate simulation-based objective functions.
L2VectorBatchManager(const Teuchos::RCP< const Teuchos::Comm< Ordinal > > &comm)
void pruneUpperActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &x, Real eps)
Set variables to zero if they correspond to the upper -active set.
std::vector< Real > x_up_
Real dx_
Definition: test_04.hpp:73
L2VectorPrimal< Real > PrimalControlVector
Definition: example_07.hpp:875
Teuchos::RCP< const std::vector< Real > > getVector() const
Definition: test_04.hpp:752
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:386
Teuchos::RCP< ROL::Vector< Real > > basis(const int i) const
Return i-th basis vector.
Definition: example_07.hpp:760
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...
void axpy(std::vector< Real > &out, const Real a, const std::vector< Real > &x, const std::vector< Real > &y) const
Real cL2_
Definition: test_04.hpp:80
Real norm() const
Returns where .
Definition: example_07.hpp:662
bool isFeasible(const ROL::Vector< Real > &x)
Check if the vector, v, is feasible.
Real dot(const ROL::Vector< Real > &x) const
Compute where .
Definition: example_07.hpp:736
void compute_residual(std::vector< Real > &r, const std::vector< Real > &u, const std::vector< Real > &z) const
Definition: example_07.hpp:322
Teuchos::RCP< std::vector< Real > > getVector()
Definition: example_07.hpp:676
Objective_BurgersControl(const Teuchos::RCP< BurgersFEM< Real > > &fem, Real x=0.0)
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:691
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:511
Teuchos::RCP< std::vector< Real > > vec_
Definition: test_04.hpp:617
void plus(const ROL::Vector< Real > &x)
Compute , where .
Definition: example_07.hpp:720
void axpy(std::vector< Real > &out, const Real a, const std::vector< Real > &x, const std::vector< Real > &y) const
Definition: example_07.hpp:90
Teuchos::RCP< ROL::Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
Definition: example_07.hpp:668
std::vector< Real > x_up_
void apply_inverse_H1(std::vector< Real > &Mu, const std::vector< Real > &u) const
Definition: example_07.hpp:293
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 ...
Teuchos::RCP< BurgersFEM< Real > > fem_
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:942
L2VectorDual(const Teuchos::RCP< std::vector< Real > > &vec, const Teuchos::RCP< BurgersFEM< Real > > &fem)
Definition: example_07.hpp:623
Teuchos::RCP< H1VectorPrimal< Real > > dual_vec_
Definition: test_04.hpp:787
Real u0_
Definition: test_04.hpp:76
void pruneActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &g, const ROL::Vector< Real > &x, Real eps)
Set variables to zero if they correspond to the -binding set.
Contains definitions of custom data types in ROL.
int dimension() const
Return dimension of the vector space.
Definition: test_04.hpp:854
void plus(const ROL::Vector< Real > &x)
Compute , where .
Definition: example_07.hpp:800
void pruneLowerActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &x, Real eps)
Set variables to zero if they correspond to the lower -active set.
void cast_vector(Teuchos::RCP< std::vector< Real > > &xvec, ROL::Vector< Real > &x) const
const std::vector< Real > getParameter(void) const
Teuchos::RCP< std::vector< Real > > vec_
Definition: test_04.hpp:704
void pruneUpperActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &g, const ROL::Vector< Real > &x, Real eps)
Set variables to zero if they correspond to the upper -binding set.
std::vector< Real >::size_type uint
Definition: example_07.hpp:881
Teuchos::RCP< std::vector< Real > > getVector()
Definition: example_07.hpp:843
Teuchos::RCP< BurgersFEM< Real > > fem_
Definition: test_04.hpp:618
Teuchos::RCP< ROL::Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
Definition: example_07.hpp:581
Real random(const Teuchos::RCP< const Teuchos::Comm< int > > &comm)
int dimension() const
Return dimension of the vector space.
Definition: test_04.hpp:600
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:358
void project(ROL::Vector< Real > &x)
Project optimization variables onto the bounds.
Teuchos::RCP< std::vector< Real > > getVector()
Definition: example_07.hpp:756
Real compute_H1_dot(const std::vector< Real > &x, const std::vector< Real > &y) const
Definition: example_07.hpp:249
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:457
virtual void zero()
Set to zero vector.
Definition: ROL_Vector.hpp:157
void scale(std::vector< Real > &u, const Real alpha=0.0) const
Definition: example_07.hpp:96
void cast_const_vector(Teuchos::RCP< const std::vector< Real > > &xvec, const ROL::Vector< Real > &x) const
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:74
Teuchos::RCP< H1VectorDual< Real > > dual_vec_
Definition: test_04.hpp:707
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 ...
H1VectorDual< Real > DualStateVector
void scale(const Real alpha)
Compute where .
Definition: example_07.hpp:562
Defines the equality constraint operator interface for simulation-based optimization.
L2VectorDual< Real > DualControlVector
Definition: example_07.hpp:876
void hessVec_22(ROL::Vector< Real > &hv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Real value(const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Compute value.
void pruneActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &x, Real eps)
Set variables to zero if they correspond to the -active set.
void plus(const ROL::Vector< Real > &x)
Compute , where .
Definition: example_07.hpp:553
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:405
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.
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:771
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:996
bool isFeasible(const ROL::Vector< Real > &x)
Check if the vector, v, is feasible.
void test_inverse_mass(std::ostream &outStream=std::cout)
Definition: example_07.hpp:212
H1BoundConstraint(std::vector< Real > &l, std::vector< Real > &u, const Teuchos::RCP< BurgersFEM< Real > > &fem, Real scale=1.0)
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:494
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:960
void apply_mass(std::vector< Real > &Mu, const std::vector< Real > &u) const
Definition: example_07.hpp:182
Teuchos::RCP< ROL::Vector< Real > > basis(const int i) const
Return i-th basis vector.
Definition: example_07.hpp:847
void apply_H1(std::vector< Real > &Mu, const std::vector< Real > &u) const
Definition: example_07.hpp:274
Real nl_
Definition: test_04.hpp:75
void pruneActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &x, Real eps)
Set variables to zero if they correspond to the -active set.
Real cH1_
Definition: test_04.hpp:79
Teuchos::RCP< BurgersFEM< Real > > fem_
Definition: test_04.hpp:785
void cast_vector(Teuchos::RCP< std::vector< Real > > &xvec, ROL::Vector< Real > &x) const
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:978
Real u1_
Definition: test_04.hpp:77
void pruneLowerActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &g, const ROL::Vector< Real > &x, Real eps)
Set variables to zero if they correspond to the lower -binding set.
void cast_const_vector(Teuchos::RCP< const std::vector< Real > > &xvec, const ROL::Vector< Real > &x) const
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:906
Real norm() const
Returns where .
Definition: example_07.hpp:742
Real compute_L2_dot(const std::vector< Real > &x, const std::vector< Real > &y) const
Definition: example_07.hpp:159
Teuchos::RCP< const std::vector< Real > > getVector() const
Definition: test_04.hpp:839
void scale(const Real alpha)
Compute where .
Definition: example_07.hpp:642
Real mesh_spacing(void) const
Definition: example_07.hpp:151
H1VectorPrimal< Real > PrimalStateVector
Definition: example_07.hpp:872
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:419
int dimension() const
Return dimension of the vector space.
Definition: test_04.hpp:687
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:604
Real norm() const
Returns where .
Definition: example_07.hpp:829
Real dot(const ROL::Vector< Real > &x) const
Compute where .
Definition: example_07.hpp:649
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:468
Real norm() const
Returns where .
Definition: example_07.hpp:575
void pruneLowerActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &x, Real eps)
Set variables to zero if they correspond to the lower -active set.
Real dot(const ROL::Vector< Real > &x) const
Compute where .
Definition: example_07.hpp:569
Teuchos::RCP< ROL::Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
Definition: example_07.hpp:835
void hessVec_21(ROL::Vector< Real > &hv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Real compute_H1_norm(const std::vector< Real > &r) const
Definition: example_07.hpp:269
std::vector< Real > x_lo_
Teuchos::RCP< L2VectorPrimal< Real > > dual_vec_
Definition: test_04.hpp:620
Teuchos::RCP< BurgersFEM< Real > > fem_
Definition: test_04.hpp:538
void hessVec_12(ROL::Vector< Real > &hv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
void setVectorToUpperBound(ROL::Vector< Real > &u)
Set the input vector to the upper bound.
H1VectorDual< Real > PrimalConstraintVector
Definition: example_07.hpp:878
Teuchos::RCP< BurgersFEM< Real > > fem_
Teuchos::RCP< const std::vector< Real > > getVector() const
Definition: test_04.hpp:672
void pruneLowerActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &g, const ROL::Vector< Real > &x, Real eps)
Set variables to zero if they correspond to the lower -binding set.
void scale(const Real alpha)
Compute where .
Definition: example_07.hpp:729
Teuchos::RCP< BurgersFEM< Real > > fem_
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.
H1VectorPrimal< Real > DualConstraintVector
Definition: example_07.hpp:879
void setVectorToLowerBound(ROL::Vector< Real > &l)
Set the input vector to the lower bound.
Teuchos::RCP< std::vector< Real > > vec_
Definition: test_04.hpp:537
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:525
Provides the interface to apply upper and lower bound constraints.
void sumAll(Real *input, Real *output, int dim)
void projection(std::vector< Real > &x)
void axpy(std::vector< Real > &out, const Real a, const std::vector< Real > &x, const std::vector< Real > &y) const
TeuchosBatchManager(const Teuchos::RCP< const Teuchos::Comm< Ordinal > > &comm)
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:858
Real nu_
Definition: test_04.hpp:74
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:924
H1VectorDual< Real > DualStateVector
Definition: example_07.hpp:873
void apply_inverse_mass(std::vector< Real > &Mu, const std::vector< Real > &u) const
Definition: example_07.hpp:199
void set_problem_data(const Real nu, const Real f, const Real u0, const Real u1)
Definition: example_07.hpp:136
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 test_inverse_H1(std::ostream &outStream=std::cout)
Definition: example_07.hpp:301
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:890
H1VectorDual(const Teuchos::RCP< std::vector< Real > > &vec, const Teuchos::RCP< BurgersFEM< Real > > &fem)
Definition: example_07.hpp:790
Real get_viscosity(void) const
Definition: example_07.hpp:143
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:102
Teuchos::RCP< BurgersFEM< Real > > fem_
Definition: test_04.hpp:881
L2VectorPrimal< Real > PrimalControlVector
void cast_vector(Teuchos::RCP< std::vector< Real > > &xvec, ROL::Vector< Real > &x) const
Teuchos::RCP< L2VectorDual< Real > > dual_vec_
Definition: test_04.hpp:540
void pruneUpperActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &g, const ROL::Vector< Real > &x, Real eps)
Set variables to zero if they correspond to the upper -binding set.
void pruneActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &g, const ROL::Vector< Real > &x, Real eps)
Set variables to zero if they correspond to the -binding set.
H1VectorPrimal< Real > PrimalStateVector
L2VectorPrimal(const Teuchos::RCP< std::vector< Real > > &vec, const Teuchos::RCP< BurgersFEM< Real > > &fem)
Definition: example_07.hpp:543
std::vector< Real > x_lo_
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:440
EqualityConstraint_BurgersControl(Teuchos::RCP< BurgersFEM< Real > > &fem, bool useHessian=true)
Definition: example_07.hpp:887
Teuchos::RCP< ROL::Vector< Real > > basis(const int i) const
Return i-th basis vector.
Definition: example_07.hpp:680
Teuchos::RCP< std::vector< Real > > vec_
Definition: test_04.hpp:784
Real f_
Definition: test_04.hpp:78
void cast_vector(Teuchos::RCP< std::vector< Real > > &xvec, ROL::Vector< Real > &x) const
std::vector< int > indices_
Real compute_L2_norm(const std::vector< Real > &r) const
Definition: example_07.hpp:177
void plus(const ROL::Vector< Real > &x)
Compute , where .
Definition: example_07.hpp:633
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:196
void sumAll(ROL::Vector< Real > &input, ROL::Vector< Real > &output)
int dimension() const
Return dimension of the vector space.
Definition: test_04.hpp:767
void project(ROL::Vector< Real > &x)
Project optimization variables onto the bounds.
void update(std::vector< Real > &u, const std::vector< Real > &s, const Real alpha=1.0) const
Definition: example_07.hpp:84
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:518
Teuchos::RCP< const std::vector< Real > > getVector() const
Definition: test_04.hpp:585
Teuchos::RCP< BurgersFEM< Real > > fem_
Definition: test_04.hpp:705
int num_dof(void) const
Definition: example_07.hpp:147
L2VectorDual< Real > DualControlVector
H1VectorPrimal(const Teuchos::RCP< std::vector< Real > > &vec, const Teuchos::RCP< BurgersFEM< Real > > &fem)
Definition: example_07.hpp:710
void sumAll(ROL::Vector< Real > &input, ROL::Vector< Real > &output)
Real dot(const ROL::Vector< Real > &x) const
Compute where .
Definition: example_07.hpp:816
Teuchos::RCP< std::vector< Real > > getVector()
Definition: example_07.hpp:589
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...
Teuchos::RCP< ROL::Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
Definition: example_07.hpp:748
L2BoundConstraint(std::vector< Real > &l, std::vector< Real > &u, const Teuchos::RCP< BurgersFEM< Real > > &fem, Real scale=1.0)
H1VectorBatchManager(const Teuchos::RCP< const Teuchos::Comm< Ordinal > > &comm)
void pruneUpperActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &x, Real eps)
Set variables to zero if they correspond to the upper -active set.
void setVectorToLowerBound(ROL::Vector< Real > &l)
Set the input vector to the lower bound.
void scale(const Real alpha)
Compute where .
Definition: example_07.hpp:809
void setVectorToUpperBound(ROL::Vector< Real > &u)
Set the input vector to the upper bound.
Teuchos::RCP< ROL::Vector< Real > > basis(const int i) const
Return i-th basis vector.
Definition: example_07.hpp:593