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