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