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"
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>
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 
1138  Real value( const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
1139  Teuchos::RCP<const std::vector<Real> > up =
1140  (Teuchos::dyn_cast<PrimalStateVector>(const_cast<ROL::Vector<Real> &>(u))).getVector();
1141 
1142 // const std::vector<Real> param
1143 // = ROL::ParametrizedObjective_SimOpt<Real>::getParameter();
1144 // fem_->set_problem_data(param[0],param[1],param[2],param[3]);
1145 // Real nu = fem_->get_viscosity();
1146 //
1147 // return 0.5*nu*fem_->compute_H1_dot(*up,*up);
1148 
1149  Real val = 0.5*((((Real)indices_[0]+1.)*(fem_->mesh_spacing())-x_)
1150  *(x_+(2.-((Real)indices_[0]+1.))*(fem_->mesh_spacing()))/(fem_->mesh_spacing())
1151  +(fem_->mesh_spacing())) * (*up)[indices_[0]];
1152  for (uint i = 1; i < indices_.size(); i++) {
1153  val += (fem_->mesh_spacing())*(*up)[indices_[i]];
1154  }
1155  return -val;
1156  }
1157 
1158  void gradient_1( ROL::Vector<Real> &g, const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
1159  Teuchos::RCP<std::vector<Real> > gp =
1160  Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<DualStateVector>(g)).getVector());
1161  Teuchos::RCP<const std::vector<Real> > up =
1162  (Teuchos::dyn_cast<PrimalStateVector>(const_cast<ROL::Vector<Real> &>(u))).getVector();
1163 
1164 // const std::vector<Real> param
1165 // = ROL::ParametrizedObjective_SimOpt<Real>::getParameter();
1166 // fem_->set_problem_data(param[0],param[1],param[2],param[3]);
1167 // Real nu = fem_->get_viscosity();
1168 //
1169 // fem_->apply_H1(*gp,*up);
1170 // g.scale(nu);
1171 
1172  g.zero();
1173  (*gp)[indices_[0]] = -0.5*((((Real)indices_[0]+1.)*(fem_->mesh_spacing())-x_)
1174  *(x_+(2.-((Real)indices_[0]+1.))*(fem_->mesh_spacing()))/(fem_->mesh_spacing())
1175  +(fem_->mesh_spacing()));
1176 
1177 
1178  for (uint i = 1; i < indices_.size(); i++) {
1179  (*gp)[indices_[i]] = -(fem_->mesh_spacing());
1180  }
1181  }
1182 
1183  void gradient_2( ROL::Vector<Real> &g, const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
1184  g.zero();
1185  }
1186 
1188  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
1189 // Teuchos::RCP<std::vector<Real> > hvp =
1190 // Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<DualStateVector>(hv)).getVector());
1191 // Teuchos::RCP<const std::vector<Real> > vp =
1192 // (Teuchos::dyn_cast<PrimalStateVector>(const_cast<ROL::Vector<Real> &>(v))).getVector();
1193 //
1194 // const std::vector<Real> param
1195 // = ROL::ParametrizedObjective_SimOpt<Real>::getParameter();
1196 // fem_->set_problem_data(param[0],param[1],param[2],param[3]);
1197 // Real nu = fem_->get_viscosity();
1198 //
1199 // fem_->apply_H1(*hvp,*vp);
1200 // hv.scale(nu);
1201 
1202  hv.zero();
1203  }
1204 
1206  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
1207  hv.zero();
1208  }
1209 
1211  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
1212  hv.zero();
1213  }
1214 
1216  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
1217  hv.zero();
1218  }
1219 };
1220 
1221 template<class Real>
1222 class L2BoundConstraint : public ROL::BoundConstraint<Real> {
1223 private:
1224  int dim_;
1225  std::vector<Real> x_lo_;
1226  std::vector<Real> x_up_;
1227  Real min_diff_;
1228  Real scale_;
1229  Teuchos::RCP<BurgersFEM<Real> > fem_;
1230 
1231  void cast_vector(Teuchos::RCP<std::vector<Real> > &xvec,
1232  ROL::Vector<Real> &x) const {
1233  try {
1234  xvec = Teuchos::rcp_const_cast<std::vector<Real> >(
1235  (Teuchos::dyn_cast<L2VectorPrimal<Real> >(x)).getVector());
1236  }
1237  catch (std::exception &e) {
1238  xvec = Teuchos::rcp_const_cast<std::vector<Real> >(
1239  (Teuchos::dyn_cast<L2VectorDual<Real> >(x)).getVector());
1240  }
1241  }
1242 
1243  void cast_const_vector(Teuchos::RCP<const std::vector<Real> > &xvec,
1244  const ROL::Vector<Real> &x) const {
1245  try {
1246  xvec = (Teuchos::dyn_cast<L2VectorPrimal<Real> >(
1247  const_cast<ROL::Vector<Real> &>(x))).getVector();
1248  }
1249  catch (std::exception &e) {
1250  xvec = (Teuchos::dyn_cast<L2VectorDual<Real> >(
1251  const_cast<ROL::Vector<Real> &>(x))).getVector();
1252  }
1253  }
1254 
1255  void axpy(std::vector<Real> &out, const Real a,
1256  const std::vector<Real> &x, const std::vector<Real> &y) const{
1257  out.resize(dim_,0.0);
1258  for (unsigned i = 0; i < dim_; i++) {
1259  out[i] = a*x[i] + y[i];
1260  }
1261  }
1262 
1263  void projection(std::vector<Real> &x) {
1264  for ( int i = 0; i < dim_; i++ ) {
1265  x[i] = std::max(x_lo_[i],std::min(x_up_[i],x[i]));
1266  }
1267  }
1268 
1269 public:
1270  L2BoundConstraint(std::vector<Real> &l, std::vector<Real> &u,
1271  const Teuchos::RCP<BurgersFEM<Real> > &fem, Real scale = 1.0)
1272  : x_lo_(l), x_up_(u), scale_(scale), fem_(fem) {
1273  dim_ = x_lo_.size();
1274  for ( int i = 0; i < dim_; i++ ) {
1275  if ( i == 0 ) {
1276  min_diff_ = x_up_[i] - x_lo_[i];
1277  }
1278  else {
1279  min_diff_ = ( (min_diff_ < (x_up_[i] - x_lo_[i])) ? min_diff_ : (x_up_[i] - x_lo_[i]) );
1280  }
1281  }
1282  min_diff_ *= 0.5;
1283  }
1284 
1285  bool isFeasible( const ROL::Vector<Real> &x ) {
1286  Teuchos::RCP<const std::vector<Real> > ex; cast_const_vector(ex,x);
1287  bool val = true;
1288  int cnt = 1;
1289  for ( int i = 0; i < dim_; i++ ) {
1290  if ( (*ex)[i] >= x_lo_[i] && (*ex)[i] <= x_up_[i] ) { cnt *= 1; }
1291  else { cnt *= 0; }
1292  }
1293  if ( cnt == 0 ) { val = false; }
1294  return val;
1295  }
1296 
1298  Teuchos::RCP<std::vector<Real> > ex; cast_vector(ex,x);
1299  projection(*ex);
1300  }
1301 
1303  Teuchos::RCP<const std::vector<Real> > ex; cast_const_vector(ex,x);
1304  Teuchos::RCP<std::vector<Real> > ev; cast_vector(ev,v);
1305  Real epsn = std::min(scale_*eps,min_diff_);
1306  for ( int i = 0; i < dim_; i++ ) {
1307  if ( ((*ex)[i] <= x_lo_[i]+epsn) ) {
1308  (*ev)[i] = 0.0;
1309  }
1310  }
1311  }
1312 
1314  Teuchos::RCP<const std::vector<Real> > ex; cast_const_vector(ex,x);
1315  Teuchos::RCP<std::vector<Real> > ev; cast_vector(ev,v);
1316  Real epsn = std::min(scale_*eps,min_diff_);
1317  for ( int i = 0; i < dim_; i++ ) {
1318  if ( ((*ex)[i] >= x_up_[i]-epsn) ) {
1319  (*ev)[i] = 0.0;
1320  }
1321  }
1322  }
1323 
1324  void pruneActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &x, Real eps) {
1325  Teuchos::RCP<const std::vector<Real> > ex; cast_const_vector(ex,x);
1326  Teuchos::RCP<std::vector<Real> > ev; cast_vector(ev,v);
1327  Real epsn = std::min(scale_*eps,min_diff_);
1328  for ( int i = 0; i < dim_; i++ ) {
1329  if ( ((*ex)[i] <= x_lo_[i]+epsn) ||
1330  ((*ex)[i] >= x_up_[i]-epsn) ) {
1331  (*ev)[i] = 0.0;
1332  }
1333  }
1334  }
1335 
1337  Teuchos::RCP<const std::vector<Real> > ex; cast_const_vector(ex,x);
1338  Teuchos::RCP<const std::vector<Real> > eg; cast_const_vector(eg,g);
1339  Teuchos::RCP<std::vector<Real> > ev; cast_vector(ev,v);
1340  Real epsn = std::min(scale_*eps,min_diff_);
1341  for ( int i = 0; i < dim_; i++ ) {
1342  if ( ((*ex)[i] <= x_lo_[i]+epsn && (*eg)[i] > 0.0) ) {
1343  (*ev)[i] = 0.0;
1344  }
1345  }
1346  }
1347 
1349  Teuchos::RCP<const std::vector<Real> > ex; cast_const_vector(ex,x);
1350  Teuchos::RCP<const std::vector<Real> > eg; cast_const_vector(eg,g);
1351  Teuchos::RCP<std::vector<Real> > ev; cast_vector(ev,v);
1352  Real epsn = std::min(scale_*eps,min_diff_);
1353  for ( int i = 0; i < dim_; i++ ) {
1354  if ( ((*ex)[i] >= x_up_[i]-epsn && (*eg)[i] < 0.0) ) {
1355  (*ev)[i] = 0.0;
1356  }
1357  }
1358  }
1359 
1360  void pruneActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &g, const ROL::Vector<Real> &x, Real eps) {
1361  Teuchos::RCP<const std::vector<Real> > ex; cast_const_vector(ex,x);
1362  Teuchos::RCP<const std::vector<Real> > eg; cast_const_vector(eg,g);
1363  Teuchos::RCP<std::vector<Real> > ev; cast_vector(ev,v);
1364  Real epsn = std::min(scale_*eps,min_diff_);
1365  for ( int i = 0; i < dim_; i++ ) {
1366  if ( ((*ex)[i] <= x_lo_[i]+epsn && (*eg)[i] > 0.0) ||
1367  ((*ex)[i] >= x_up_[i]-epsn && (*eg)[i] < 0.0) ) {
1368  (*ev)[i] = 0.0;
1369  }
1370  }
1371  }
1372 
1374  Teuchos::RCP<std::vector<Real> > us = Teuchos::rcp( new std::vector<Real>(dim_,0.0) );
1375  us->assign(x_up_.begin(),x_up_.end());
1376  Teuchos::RCP<ROL::Vector<Real> > up = Teuchos::rcp( new L2VectorPrimal<Real>(us,fem_) );
1377  u.set(*up);
1378  }
1379 
1381  Teuchos::RCP<std::vector<Real> > ls = Teuchos::rcp( new std::vector<Real>(dim_,0.0) );
1382  ls->assign(x_lo_.begin(),x_lo_.end());
1383  Teuchos::RCP<ROL::Vector<Real> > lp = Teuchos::rcp( new L2VectorPrimal<Real>(ls,fem_) );
1384  l.set(*lp);
1385  }
1386 };
1387 
1388 template<class Real>
1389 class H1BoundConstraint : public ROL::BoundConstraint<Real> {
1390 private:
1391  int dim_;
1392  std::vector<Real> x_lo_;
1393  std::vector<Real> x_up_;
1394  Real min_diff_;
1395  Real scale_;
1396  Teuchos::RCP<BurgersFEM<Real> > fem_;
1397 
1398  void cast_vector(Teuchos::RCP<std::vector<Real> > &xvec,
1399  ROL::Vector<Real> &x) const {
1400  try {
1401  xvec = Teuchos::rcp_const_cast<std::vector<Real> >(
1402  (Teuchos::dyn_cast<H1VectorPrimal<Real> >(x)).getVector());
1403  }
1404  catch (std::exception &e) {
1405  xvec = Teuchos::rcp_const_cast<std::vector<Real> >(
1406  (Teuchos::dyn_cast<H1VectorDual<Real> >(x)).getVector());
1407  }
1408  }
1409 
1410  void cast_const_vector(Teuchos::RCP<const std::vector<Real> > &xvec,
1411  const ROL::Vector<Real> &x) const {
1412  try {
1413  xvec = (Teuchos::dyn_cast<H1VectorPrimal<Real> >(
1414  const_cast<ROL::Vector<Real> &>(x))).getVector();
1415  }
1416  catch (std::exception &e) {
1417  xvec = (Teuchos::dyn_cast<H1VectorDual<Real> >(
1418  const_cast<ROL::Vector<Real> &>(x))).getVector();
1419  }
1420  }
1421 
1422  void axpy(std::vector<Real> &out, const Real a,
1423  const std::vector<Real> &x, const std::vector<Real> &y) const{
1424  out.resize(dim_,0.0);
1425  for (unsigned i = 0; i < dim_; i++) {
1426  out[i] = a*x[i] + y[i];
1427  }
1428  }
1429 
1430  void projection(std::vector<Real> &x) {
1431  for ( int i = 0; i < dim_; i++ ) {
1432  x[i] = std::max(x_lo_[i],std::min(x_up_[i],x[i]));
1433  }
1434  }
1435 
1436 public:
1437  H1BoundConstraint(std::vector<Real> &l, std::vector<Real> &u,
1438  const Teuchos::RCP<BurgersFEM<Real> > &fem, Real scale = 1.0)
1439  : x_lo_(l), x_up_(u), scale_(scale), fem_(fem) {
1440  dim_ = x_lo_.size();
1441  for ( int i = 0; i < dim_; i++ ) {
1442  if ( i == 0 ) {
1443  min_diff_ = x_up_[i] - x_lo_[i];
1444  }
1445  else {
1446  min_diff_ = ( (min_diff_ < (x_up_[i] - x_lo_[i])) ? min_diff_ : (x_up_[i] - x_lo_[i]) );
1447  }
1448  }
1449  min_diff_ *= 0.5;
1450  }
1451 
1452  bool isFeasible( const ROL::Vector<Real> &x ) {
1453  Teuchos::RCP<const std::vector<Real> > ex; cast_const_vector(ex,x);
1454  bool val = true;
1455  int cnt = 1;
1456  for ( int i = 0; i < dim_; i++ ) {
1457  if ( (*ex)[i] >= x_lo_[i] && (*ex)[i] <= x_up_[i] ) { cnt *= 1; }
1458  else { cnt *= 0; }
1459  }
1460  if ( cnt == 0 ) { val = false; }
1461  return val;
1462  }
1463 
1465  Teuchos::RCP<std::vector<Real> > ex; cast_vector(ex,x);
1466  projection(*ex);
1467  }
1468 
1470  Teuchos::RCP<const std::vector<Real> > ex; cast_const_vector(ex,x);
1471  Teuchos::RCP<std::vector<Real> > ev; cast_vector(ev,v);
1472  Real epsn = std::min(scale_*eps,min_diff_);
1473  for ( int i = 0; i < dim_; i++ ) {
1474  if ( ((*ex)[i] <= x_lo_[i]+epsn) ) {
1475  (*ev)[i] = 0.0;
1476  }
1477  }
1478  }
1479 
1481  Teuchos::RCP<const std::vector<Real> > ex; cast_const_vector(ex,x);
1482  Teuchos::RCP<std::vector<Real> > ev; cast_vector(ev,v);
1483  Real epsn = std::min(scale_*eps,min_diff_);
1484  for ( int i = 0; i < dim_; i++ ) {
1485  if ( ((*ex)[i] >= x_up_[i]-epsn) ) {
1486  (*ev)[i] = 0.0;
1487  }
1488  }
1489  }
1490 
1491  void pruneActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &x, Real eps) {
1492  Teuchos::RCP<const std::vector<Real> > ex; cast_const_vector(ex,x);
1493  Teuchos::RCP<std::vector<Real> > ev; cast_vector(ev,v);
1494  Real epsn = std::min(scale_*eps,min_diff_);
1495  for ( int i = 0; i < dim_; i++ ) {
1496  if ( ((*ex)[i] <= x_lo_[i]+epsn) ||
1497  ((*ex)[i] >= x_up_[i]-epsn) ) {
1498  (*ev)[i] = 0.0;
1499  }
1500  }
1501  }
1502 
1504  Teuchos::RCP<const std::vector<Real> > ex; cast_const_vector(ex,x);
1505  Teuchos::RCP<const std::vector<Real> > eg; cast_const_vector(eg,g);
1506  Teuchos::RCP<std::vector<Real> > ev; cast_vector(ev,v);
1507  Real epsn = std::min(scale_*eps,min_diff_);
1508  for ( int i = 0; i < dim_; i++ ) {
1509  if ( ((*ex)[i] <= x_lo_[i]+epsn && (*eg)[i] > 0.0) ) {
1510  (*ev)[i] = 0.0;
1511  }
1512  }
1513  }
1514 
1516  Teuchos::RCP<const std::vector<Real> > ex; cast_const_vector(ex,x);
1517  Teuchos::RCP<const std::vector<Real> > eg; cast_const_vector(eg,g);
1518  Teuchos::RCP<std::vector<Real> > ev; cast_vector(ev,v);
1519  Real epsn = std::min(scale_*eps,min_diff_);
1520  for ( int i = 0; i < dim_; i++ ) {
1521  if ( ((*ex)[i] >= x_up_[i]-epsn && (*eg)[i] < 0.0) ) {
1522  (*ev)[i] = 0.0;
1523  }
1524  }
1525  }
1526 
1527  void pruneActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &g, const ROL::Vector<Real> &x, Real eps) {
1528  Teuchos::RCP<const std::vector<Real> > ex; cast_const_vector(ex,x);
1529  Teuchos::RCP<const std::vector<Real> > eg; cast_const_vector(eg,g);
1530  Teuchos::RCP<std::vector<Real> > ev; cast_vector(ev,v);
1531  Real epsn = std::min(scale_*eps,min_diff_);
1532  for ( int i = 0; i < dim_; i++ ) {
1533  if ( ((*ex)[i] <= x_lo_[i]+epsn && (*eg)[i] > 0.0) ||
1534  ((*ex)[i] >= x_up_[i]-epsn && (*eg)[i] < 0.0) ) {
1535  (*ev)[i] = 0.0;
1536  }
1537  }
1538  }
1539 
1541  Teuchos::RCP<std::vector<Real> > us = Teuchos::rcp( new std::vector<Real>(dim_,0.0) );
1542  us->assign(x_up_.begin(),x_up_.end());
1543  Teuchos::RCP<ROL::Vector<Real> > up = Teuchos::rcp( new H1VectorPrimal<Real>(us,fem_) );
1544  u.set(*up);
1545  }
1546 
1548  Teuchos::RCP<std::vector<Real> > ls = Teuchos::rcp( new std::vector<Real>(dim_,0.0) );
1549  ls->assign(x_lo_.begin(),x_lo_.end());
1550  Teuchos::RCP<ROL::Vector<Real> > lp = Teuchos::rcp( new H1VectorPrimal<Real>(ls,fem_) );
1551  l.set(*lp);
1552  }
1553 };
1554 
1555 template<class Real, class Ordinal>
1556 class L2VectorBatchManager : public ROL::TeuchosBatchManager<Real,Ordinal> {
1557 private:
1558  void cast_vector(Teuchos::RCP<std::vector<Real> > &xvec,
1559  ROL::Vector<Real> &x) const {
1560  try {
1561  xvec = Teuchos::rcp_const_cast<std::vector<Real> >(
1562  (Teuchos::dyn_cast<L2VectorPrimal<Real> >(x)).getVector());
1563  }
1564  catch (std::exception &e) {
1565  xvec = Teuchos::rcp_const_cast<std::vector<Real> >(
1566  (Teuchos::dyn_cast<L2VectorDual<Real> >(x)).getVector());
1567  }
1568  }
1569 
1570 public:
1571  L2VectorBatchManager(const Teuchos::RCP<const Teuchos::Comm<Ordinal> > &comm)
1572  : ROL::TeuchosBatchManager<Real,Ordinal>(comm) {}
1574  Teuchos::RCP<std::vector<Real> > input_ptr;
1575  cast_vector(input_ptr,input);
1576  int dim_i = input_ptr->size();
1577  Teuchos::RCP<std::vector<Real> > output_ptr;
1578  cast_vector(output_ptr,output);
1579  int dim_o = output_ptr->size();
1580  if ( dim_i != dim_o ) {
1581  std::cout << "L2VectorBatchManager: DIMENSION MISMATCH ON RANK "
1583  }
1584  else {
1585  ROL::TeuchosBatchManager<Real,Ordinal>::sumAll(&(*input_ptr)[0],&(*output_ptr)[0],dim_i);
1586  }
1587  }
1588 };
1589 
1590 template<class Real, class Ordinal>
1591 class H1VectorBatchManager : public ROL::TeuchosBatchManager<Real,Ordinal> {
1592 private:
1593  void cast_vector(Teuchos::RCP<std::vector<Real> > &xvec,
1594  ROL::Vector<Real> &x) const {
1595  try {
1596  xvec = Teuchos::rcp_const_cast<std::vector<Real> >(
1597  (Teuchos::dyn_cast<H1VectorPrimal<Real> >(x)).getVector());
1598  }
1599  catch (std::exception &e) {
1600  xvec = Teuchos::rcp_const_cast<std::vector<Real> >(
1601  (Teuchos::dyn_cast<H1VectorDual<Real> >(x)).getVector());
1602  }
1603  }
1604 
1605 public:
1606  H1VectorBatchManager(const Teuchos::RCP<const Teuchos::Comm<Ordinal> > &comm)
1607  : ROL::TeuchosBatchManager<Real,Ordinal>(comm) {}
1609  Teuchos::RCP<std::vector<Real> > input_ptr;
1610  cast_vector(input_ptr,input);
1611  int dim_i = input_ptr->size();
1612  Teuchos::RCP<std::vector<Real> > output_ptr;
1613  cast_vector(output_ptr,output);
1614  int dim_o = output_ptr->size();
1615  if ( dim_i != dim_o ) {
1616  std::cout << "H1VectorBatchManager: DIMENSION MISMATCH ON RANK "
1618  }
1619  else {
1620  ROL::TeuchosBatchManager<Real,Ordinal>::sumAll(&(*input_ptr)[0],&(*output_ptr)[0],dim_i);
1621  }
1622  }
1623 };
1624 
1625 template<class Real>
1626 Real random(const Teuchos::RCP<const Teuchos::Comm<int> > &comm) {
1627  Real val = 0.0;
1628  if ( Teuchos::rank<int>(*comm)==0 ) {
1629  val = (Real)rand()/(Real)RAND_MAX;
1630  }
1631  Teuchos::broadcast<int,Real>(*comm,0,1,&val);
1632  return val;
1633 }
BurgersFEM(int nx=128, Real nl=1.0, Real cH1=1.0, Real cL2=1.0)
Definition: example_07.hpp:133
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.
void cast_vector(Teuchos::RCP< std::vector< Real > > &xvec, ROL::Vector< Real > &x) const
Real dx_
Definition: test_04.hpp:73
L2VectorPrimal< Real > PrimalControlVector
Definition: example_07.hpp:875
Real compute_H1_dot(const std::vector< Real > &x, const std::vector< Real > &y) const
Definition: example_07.hpp:249
void applyAdjointHessian_21(ROL::Vector< Real > &ahwv, const ROL::Vector< Real > &w, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Apply the adjoint of the partial constraint Hessian at , , to vector in direction ...
Real norm() const
Returns where .
Definition: example_07.hpp:742
Real cL2_
Definition: test_04.hpp:80
void cast_const_vector(Teuchos::RCP< const std::vector< Real > > &xvec, const ROL::Vector< Real > &x) const
void apply_mass(std::vector< Real > &Mu, const std::vector< Real > &u) const
Definition: example_07.hpp:182
Teuchos::RCP< const std::vector< Real > > getVector() const
Definition: test_04.hpp:672
bool isFeasible(const ROL::Vector< Real > &x)
Check if the vector, v, is feasible.
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< std::vector< Real > > getVector()
Definition: example_07.hpp:676
Teuchos::RCP< ROL::Vector< Real > > basis(const int i) const
Return i-th basis vector.
Definition: example_07.hpp:593
Objective_BurgersControl(const Teuchos::RCP< BurgersFEM< Real > > &fem, Real x=0.0)
Real dot(const ROL::Vector< Real > &x) const
Compute where .
Definition: example_07.hpp:649
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< const std::vector< Real > > getVector() const
Definition: test_04.hpp:585
int dimension() const
Return dimension of the vector space.
Definition: example_07.hpp:600
Real get_viscosity(void) const
Definition: example_07.hpp:143
void plus(const ROL::Vector< Real > &x)
Compute , where .
Definition: example_07.hpp:720
void applyAdjointHessian_22(ROL::Vector< Real > &ahwv, const ROL::Vector< Real > &w, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Apply the adjoint of the partial constraint Hessian at , , to vector in direction ...
void applyInverseJacobian_1(ROL::Vector< Real > &ijv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Apply the inverse partial constraint Jacobian at , , to the vector .
Definition: example_07.hpp:942
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
L2VectorDual(const Teuchos::RCP< std::vector< Real > > &vec, const Teuchos::RCP< BurgersFEM< Real > > &fem)
Definition: example_07.hpp:623
Teuchos::RCP< ROL::Vector< Real > > basis(const int i) const
Return i-th basis vector.
Definition: example_07.hpp:847
Real u0_
Definition: test_04.hpp:76
Real norm() const
Returns where .
Definition: example_07.hpp:575
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.
Real dot(const ROL::Vector< Real > &x) const
Compute where .
Definition: example_07.hpp:816
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.
Teuchos::RCP< const std::vector< Real > > getVector() const
Definition: test_04.hpp:752
Real compute_H1_norm(const std::vector< Real > &r) const
Definition: example_07.hpp:269
Real norm() const
Returns where .
Definition: example_07.hpp:662
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
Real random(const Teuchos::RCP< const Teuchos::Comm< int > > &comm)
void project(ROL::Vector< Real > &x)
Project optimization variables onto the bounds.
Teuchos::RCP< std::vector< Real > > getVector()
Definition: example_07.hpp:756
virtual void zero()
Set to zero vector.
Definition: ROL_Vector.hpp:157
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:74
void applyAdjointHessian_11(ROL::Vector< Real > &ahwv, const ROL::Vector< Real > &w, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Apply the adjoint of the partial constraint Hessian at , , to vector in direction ...
H1VectorDual< Real > DualStateVector
void scale(const Real alpha)
Compute where .
Definition: example_07.hpp:562
L2VectorDual< Real > DualControlVector
Definition: example_07.hpp:876
Real dot(const ROL::Vector< Real > &x) const
Compute where .
Definition: example_07.hpp:736
int num_dof(void) const
Definition: example_07.hpp:147
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 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:151
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.
int dimension() const
Return dimension of the vector space.
Definition: example_07.hpp:767
void test_inverse_mass(std::ostream &outStream=std::cout)
Definition: example_07.hpp:212
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
H1BoundConstraint(std::vector< Real > &l, std::vector< Real > &u, const Teuchos::RCP< BurgersFEM< Real > > &fem, Real scale=1.0)
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
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
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 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
Teuchos::RCP< ROL::Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
Definition: example_07.hpp:748
void scale(const Real alpha)
Compute where .
Definition: example_07.hpp:642
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
H1VectorPrimal< Real > PrimalStateVector
Definition: example_07.hpp:872
Teuchos::RCP< ROL::Vector< Real > > basis(const int i) const
Return i-th basis vector.
Definition: example_07.hpp:760
Real compute_L2_dot(const std::vector< Real > &x, const std::vector< Real > &y) const
Definition: example_07.hpp:159
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 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:102
void cast_vector(Teuchos::RCP< std::vector< Real > > &xvec, ROL::Vector< Real > &x) const
void apply_inverse_mass(std::vector< Real > &Mu, const std::vector< Real > &u) const
Definition: example_07.hpp:199
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)
void setVectorToUpperBound(ROL::Vector< Real > &u)
Set the input vector to the upper bound.
H1VectorDual< Real > PrimalConstraintVector
Definition: example_07.hpp:878
void update(std::vector< Real > &u, const std::vector< Real > &s, const Real alpha=1.0) const
Definition: example_07.hpp:84
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
void projection(std::vector< Real > &x)
void gradient_2(ROL::Vector< Real > &g, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Compute gradient with respect to second component.
void cast_vector(Teuchos::RCP< std::vector< Real > > &xvec, ROL::Vector< Real > &x) const
H1VectorPrimal< Real > DualConstraintVector
Definition: example_07.hpp:879
void setVectorToLowerBound(ROL::Vector< Real > &l)
Set the input vector to the lower bound.
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 compute_residual(std::vector< Real > &r, const std::vector< Real > &u, const std::vector< Real > &z) const
Definition: example_07.hpp:322
void projection(std::vector< Real > &x)
Real nu_
Definition: test_04.hpp:74
void cast_const_vector(Teuchos::RCP< const std::vector< Real > > &xvec, const ROL::Vector< Real > &x) const
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 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 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
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
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
H1VectorDual(const Teuchos::RCP< std::vector< Real > > &vec, const Teuchos::RCP< BurgersFEM< Real > > &fem)
Definition: example_07.hpp:790
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 dot(const ROL::Vector< Real > &x) const
Compute where .
Definition: example_07.hpp:569
L2VectorPrimal< Real > PrimalControlVector
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
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 axpy(std::vector< Real > &out, const Real a, const std::vector< Real > &x, const std::vector< Real > &y) const
void cast_vector(Teuchos::RCP< std::vector< Real > > &xvec, ROL::Vector< Real > &x) const
EqualityConstraint_BurgersControl(Teuchos::RCP< BurgersFEM< Real > > &fem, bool useHessian=true)
Definition: example_07.hpp:887
Teuchos::RCP< ROL::Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
Definition: example_07.hpp:581
Real f_
Definition: test_04.hpp:78
std::vector< int > indices_
void apply_H1(std::vector< Real > &Mu, const std::vector< Real > &u) const
Definition: example_07.hpp:274
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)
void apply_inverse_H1(std::vector< Real > &Mu, const std::vector< Real > &u) const
Definition: example_07.hpp:293
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:419
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
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
int dimension() const
Return dimension of the vector space.
Definition: example_07.hpp:854
L2VectorDual< Real > DualControlVector
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
Teuchos::RCP< ROL::Vector< Real > > basis(const int i) const
Return i-th basis vector.
Definition: example_07.hpp:680
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 norm() const
Returns where .
Definition: example_07.hpp:829
void scale(std::vector< Real > &u, const Real alpha=0.0) const
Definition: example_07.hpp:96
Teuchos::RCP< std::vector< Real > > getVector()
Definition: example_07.hpp:589
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
Teuchos::RCP< ROL::Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
Definition: example_07.hpp:835
void applyAdjointHessian_12(ROL::Vector< Real > &ahwv, const ROL::Vector< Real > &w, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Apply the adjoint of the partial constraint Hessian at , , to vector in direction ...
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.
Teuchos::RCP< const std::vector< Real > > getVector() const
Definition: test_04.hpp:839
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
int dimension() const
Return dimension of the vector space.
Definition: example_07.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 compute_L2_norm(const std::vector< Real > &r) const
Definition: example_07.hpp:177
Teuchos::RCP< ROL::Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
Definition: example_07.hpp:668
void setVectorToUpperBound(ROL::Vector< Real > &u)
Set the input vector to the upper bound.