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"
52 #include "ROL_Constraint_SimOpt.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  ROL::Ptr<std::vector<Real> > vec_;
530  ROL::Ptr<BurgersFEM<Real> > fem_;
531 
532  mutable ROL::Ptr<L2VectorDual<Real> > dual_vec_;
533 
534 public:
535  L2VectorPrimal(const ROL::Ptr<std::vector<Real> > & vec,
536  const ROL::Ptr<BurgersFEM<Real> > &fem)
537  : vec_(vec), fem_(fem), dual_vec_(ROL::nullPtr) {}
538 
539  void set( const ROL::Vector<Real> &x ) {
540  const L2VectorPrimal &ex = dynamic_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 = dynamic_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 = dynamic_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  ROL::Ptr<ROL::Vector<Real> > clone() const {
574  return ROL::makePtr<L2VectorPrimal>( ROL::makePtr<std::vector<Real>>(vec_->size(),0.0),fem_);
575  }
576 
577  ROL::Ptr<const std::vector<Real> > getVector() const {
578  return vec_;
579  }
580 
581  ROL::Ptr<std::vector<Real> > getVector() {
582  return vec_;
583  }
584 
585  ROL::Ptr<ROL::Vector<Real> > basis( const int i ) const {
586  ROL::Ptr<L2VectorPrimal> e
587  = ROL::makePtr<L2VectorPrimal>( ROL::makePtr<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_ = ROL::makePtr<L2VectorDual<Real>>(
598  ROL::makePtr<std::vector<Real>>(*vec_),fem_);
599 
600  fem_->apply_mass(*(ROL::constPtrCast<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  ROL::Ptr<std::vector<Real> > vec_;
610  ROL::Ptr<BurgersFEM<Real> > fem_;
611 
612  mutable ROL::Ptr<L2VectorPrimal<Real> > dual_vec_;
613 
614 public:
615  L2VectorDual(const ROL::Ptr<std::vector<Real> > & vec,
616  const ROL::Ptr<BurgersFEM<Real> > &fem)
617  : vec_(vec), fem_(fem), dual_vec_(ROL::nullPtr) {}
618 
619  void set( const ROL::Vector<Real> &x ) {
620  const L2VectorDual &ex = dynamic_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 = dynamic_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 = dynamic_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  ROL::Ptr<ROL::Vector<Real> > clone() const {
661  return ROL::makePtr<L2VectorDual>( ROL::makePtr<std::vector<Real>>(vec_->size(),0.0),fem_);
662  }
663 
664  ROL::Ptr<const std::vector<Real> > getVector() const {
665  return vec_;
666  }
667 
668  ROL::Ptr<std::vector<Real> > getVector() {
669  return vec_;
670  }
671 
672  ROL::Ptr<ROL::Vector<Real> > basis( const int i ) const {
673  ROL::Ptr<L2VectorDual> e
674  = ROL::makePtr<L2VectorDual>( ROL::makePtr<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_ = ROL::makePtr<L2VectorPrimal<Real>>(
685  ROL::makePtr<std::vector<Real>>(*vec_),fem_);
686 
687  fem_->apply_inverse_mass(*(ROL::constPtrCast<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  ROL::Ptr<std::vector<Real> > vec_;
697  ROL::Ptr<BurgersFEM<Real> > fem_;
698 
699  mutable ROL::Ptr<H1VectorDual<Real> > dual_vec_;
700 
701 public:
702  H1VectorPrimal(const ROL::Ptr<std::vector<Real> > & vec,
703  const ROL::Ptr<BurgersFEM<Real> > &fem)
704  : vec_(vec), fem_(fem), dual_vec_(ROL::nullPtr) {}
705 
706  void set( const ROL::Vector<Real> &x ) {
707  const H1VectorPrimal &ex = dynamic_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 = dynamic_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 = dynamic_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  ROL::Ptr<ROL::Vector<Real> > clone() const {
741  return ROL::makePtr<H1VectorPrimal>( ROL::makePtr<std::vector<Real>>(vec_->size(),0.0),fem_);
742  }
743 
744  ROL::Ptr<const std::vector<Real> > getVector() const {
745  return vec_;
746  }
747 
748  ROL::Ptr<std::vector<Real> > getVector() {
749  return vec_;
750  }
751 
752  ROL::Ptr<ROL::Vector<Real> > basis( const int i ) const {
753  ROL::Ptr<H1VectorPrimal> e
754  = ROL::makePtr<H1VectorPrimal>( ROL::makePtr<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_ = ROL::makePtr<H1VectorDual<Real>>(
765  ROL::makePtr<std::vector<Real>>(*vec_),fem_);
766 
767  fem_->apply_H1(*(ROL::constPtrCast<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  ROL::Ptr<std::vector<Real> > vec_;
777  ROL::Ptr<BurgersFEM<Real> > fem_;
778 
779  mutable ROL::Ptr<H1VectorPrimal<Real> > dual_vec_;
780 
781 public:
782  H1VectorDual(const ROL::Ptr<std::vector<Real> > & vec,
783  const ROL::Ptr<BurgersFEM<Real> > &fem)
784  : vec_(vec), fem_(fem), dual_vec_(ROL::nullPtr) {}
785 
786  void set( const ROL::Vector<Real> &x ) {
787  const H1VectorDual &ex = dynamic_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 = dynamic_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 = dynamic_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  ROL::Ptr<ROL::Vector<Real> > clone() const {
828  return ROL::makePtr<H1VectorDual>( ROL::makePtr<std::vector<Real>>(vec_->size(),0.0),fem_);
829  }
830 
831  ROL::Ptr<const std::vector<Real> > getVector() const {
832  return vec_;
833  }
834 
835  ROL::Ptr<std::vector<Real> > getVector() {
836  return vec_;
837  }
838 
839  ROL::Ptr<ROL::Vector<Real> > basis( const int i ) const {
840  ROL::Ptr<H1VectorDual> e
841  = ROL::makePtr<H1VectorDual>( ROL::makePtr<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_ = ROL::makePtr<H1VectorPrimal<Real>>(
852  ROL::makePtr<std::vector<Real>>(*vec_),fem_);
853 
854  fem_->apply_inverse_H1(*(ROL::constPtrCast<std::vector<Real> >(dual_vec_->getVector())),*vec_);
855  return *dual_vec_;
856  }
857 
858 };
859 
860 template<class Real>
861 class Constraint_BurgersControl : public ROL::Constraint_SimOpt<Real> {
862 private:
863 
864  typedef H1VectorPrimal<Real> PrimalStateVector;
865  typedef H1VectorDual<Real> DualStateVector;
866 
867  typedef L2VectorPrimal<Real> PrimalControlVector;
868  typedef L2VectorDual<Real> DualControlVector;
869 
870  typedef H1VectorDual<Real> PrimalConstraintVector;
871  typedef H1VectorPrimal<Real> DualConstraintVector;
872 
873  ROL::Ptr<BurgersFEM<Real> > fem_;
874  bool useHessian_;
875 
876 public:
877  Constraint_BurgersControl(ROL::Ptr<BurgersFEM<Real> > &fem, bool useHessian = true)
878  : fem_(fem), useHessian_(useHessian) {}
879 
880  void value(ROL::Vector<Real> &c, const ROL::Vector<Real> &u,
881  const ROL::Vector<Real> &z, Real &tol) {
882  ROL::Ptr<std::vector<Real> > cp =
883  dynamic_cast<PrimalConstraintVector&>(c).getVector();
884  ROL::Ptr<const std::vector<Real> > up =
885  dynamic_cast<const PrimalStateVector&>(u).getVector();
886  ROL::Ptr<const std::vector<Real> > zp =
887  dynamic_cast<const PrimalControlVector&>(z).getVector();
888  fem_->compute_residual(*cp,*up,*zp);
889  }
890 
891  using ROL::Constraint_SimOpt<Real>::value;
892 
893  void applyJacobian_1(ROL::Vector<Real> &jv, const ROL::Vector<Real> &v, const ROL::Vector<Real> &u,
894  const ROL::Vector<Real> &z, Real &tol) {
895  ROL::Ptr<std::vector<Real> > jvp =
896  dynamic_cast<PrimalConstraintVector&>(jv).getVector();
897  ROL::Ptr<const std::vector<Real> > vp =
898  dynamic_cast<const PrimalStateVector&>(v).getVector();
899  ROL::Ptr<const std::vector<Real> > up =
900  dynamic_cast<const PrimalStateVector&>(u).getVector();
901  ROL::Ptr<const std::vector<Real> > zp =
902  dynamic_cast<const PrimalControlVector&>(z).getVector();
903  fem_->apply_pde_jacobian(*jvp,*vp,*up,*zp);
904  }
905 
906  void applyJacobian_2(ROL::Vector<Real> &jv, const ROL::Vector<Real> &v, const ROL::Vector<Real> &u,
907  const ROL::Vector<Real> &z, Real &tol) {
908  ROL::Ptr<std::vector<Real> > jvp =
909  dynamic_cast<PrimalConstraintVector&>(jv).getVector();
910  ROL::Ptr<const std::vector<Real> > vp =
911  dynamic_cast<const PrimalControlVector&>(v).getVector();
912  ROL::Ptr<const std::vector<Real> > up =
913  dynamic_cast<const PrimalStateVector&>(u).getVector();
914  ROL::Ptr<const std::vector<Real> > zp =
915  dynamic_cast<const PrimalControlVector&>(z).getVector();
916  fem_->apply_control_jacobian(*jvp,*vp,*up,*zp);
917  }
918 
919  void applyInverseJacobian_1(ROL::Vector<Real> &ijv, const ROL::Vector<Real> &v, const ROL::Vector<Real> &u,
920  const ROL::Vector<Real> &z, Real &tol) {
921  ROL::Ptr<std::vector<Real> > ijvp =
922  dynamic_cast<PrimalStateVector&>(ijv).getVector();
923  ROL::Ptr<const std::vector<Real> > vp =
924  dynamic_cast<const PrimalConstraintVector&>(v).getVector();
925  ROL::Ptr<const std::vector<Real> > up =
926  dynamic_cast<const PrimalStateVector&>(u).getVector();
927  ROL::Ptr<const std::vector<Real> > zp =
928  dynamic_cast<const PrimalControlVector&>(z).getVector();
929  fem_->apply_inverse_pde_jacobian(*ijvp,*vp,*up,*zp);
930  }
931 
932  void applyAdjointJacobian_1(ROL::Vector<Real> &ajv, const ROL::Vector<Real> &v, const ROL::Vector<Real> &u,
933  const ROL::Vector<Real> &z, Real &tol) {
934  ROL::Ptr<std::vector<Real> > jvp =
935  dynamic_cast<DualStateVector&>(ajv).getVector();
936  ROL::Ptr<const std::vector<Real> > vp =
937  dynamic_cast<const DualConstraintVector&>(v).getVector();
938  ROL::Ptr<const std::vector<Real> > up =
939  dynamic_cast<const PrimalStateVector&>(u).getVector();
940  ROL::Ptr<const std::vector<Real> > zp =
941  dynamic_cast<const PrimalControlVector&>(z).getVector();
942  fem_->apply_adjoint_pde_jacobian(*jvp,*vp,*up,*zp);
943  }
944 
945  void applyAdjointJacobian_2(ROL::Vector<Real> &jv, const ROL::Vector<Real> &v, const ROL::Vector<Real> &u,
946  const ROL::Vector<Real> &z, Real &tol) {
947  ROL::Ptr<std::vector<Real> > jvp =
948  dynamic_cast<DualControlVector&>(jv).getVector();
949  ROL::Ptr<const std::vector<Real> > vp =
950  dynamic_cast<const DualConstraintVector&>(v).getVector();
951  ROL::Ptr<const std::vector<Real> > up =
952  dynamic_cast<const PrimalStateVector&>(u).getVector();
953  ROL::Ptr<const std::vector<Real> > zp =
954  dynamic_cast<const PrimalControlVector&>(z).getVector();
955  fem_->apply_adjoint_control_jacobian(*jvp,*vp,*up,*zp);
956  }
957 
958  void applyInverseAdjointJacobian_1(ROL::Vector<Real> &iajv, const ROL::Vector<Real> &v,
959  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
960  ROL::Ptr<std::vector<Real> > iajvp =
961  dynamic_cast<DualConstraintVector&>(iajv).getVector();
962  ROL::Ptr<const std::vector<Real> > vp =
963  dynamic_cast<const DualStateVector&>(v).getVector();
964  ROL::Ptr<const std::vector<Real> > up =
965  dynamic_cast<const PrimalStateVector&>(u).getVector();
966  ROL::Ptr<const std::vector<Real> > zp =
967  dynamic_cast<const PrimalControlVector&>(z).getVector();
968  fem_->apply_inverse_adjoint_pde_jacobian(*iajvp,*vp,*up,*zp);
969  }
970 
971  void applyAdjointHessian_11(ROL::Vector<Real> &ahwv, const ROL::Vector<Real> &w, const ROL::Vector<Real> &v,
972  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
973  if ( useHessian_ ) {
974  ROL::Ptr<std::vector<Real> > ahwvp =
975  dynamic_cast<DualStateVector&>(ahwv).getVector();
976  ROL::Ptr<const std::vector<Real> > wp =
977  dynamic_cast<const DualConstraintVector&>(w).getVector();
978  ROL::Ptr<const std::vector<Real> > vp =
979  dynamic_cast<const PrimalStateVector&>(v).getVector();
980  ROL::Ptr<const std::vector<Real> > up =
981  dynamic_cast<const PrimalStateVector&>(u).getVector();
982  ROL::Ptr<const std::vector<Real> > zp =
983  dynamic_cast<const PrimalControlVector&>(z).getVector();
984  fem_->apply_adjoint_pde_hessian(*ahwvp,*wp,*vp,*up,*zp);
985  }
986  else {
987  ahwv.zero();
988  }
989  }
990 
991  void applyAdjointHessian_12(ROL::Vector<Real> &ahwv, const ROL::Vector<Real> &w, const ROL::Vector<Real> &v,
992  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
993  if ( useHessian_ ) {
994  ROL::Ptr<std::vector<Real> > ahwvp =
995  dynamic_cast<DualControlVector&>(ahwv).getVector();
996  ROL::Ptr<const std::vector<Real> > wp =
997  dynamic_cast<const DualConstraintVector&>(w).getVector();
998  ROL::Ptr<const std::vector<Real> > vp =
999  dynamic_cast<const PrimalStateVector&>(v).getVector();
1000  ROL::Ptr<const std::vector<Real> > up =
1001  dynamic_cast<const PrimalStateVector&>(u).getVector();
1002  ROL::Ptr<const std::vector<Real> > zp =
1003  dynamic_cast<const PrimalControlVector&>(z).getVector();
1004  fem_->apply_adjoint_control_pde_hessian(*ahwvp,*wp,*vp,*up,*zp);
1005  }
1006  else {
1007  ahwv.zero();
1008  }
1009  }
1010  void applyAdjointHessian_21(ROL::Vector<Real> &ahwv, const ROL::Vector<Real> &w, const ROL::Vector<Real> &v,
1011  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
1012  if ( useHessian_ ) {
1013  ROL::Ptr<std::vector<Real> > ahwvp =
1014  dynamic_cast<DualStateVector&>(ahwv).getVector();
1015  ROL::Ptr<const std::vector<Real> > wp =
1016  dynamic_cast<const DualConstraintVector&>(w).getVector();
1017  ROL::Ptr<const std::vector<Real> > vp =
1018  dynamic_cast<const PrimalControlVector&>(v).getVector();
1019  ROL::Ptr<const std::vector<Real> > up =
1020  dynamic_cast<const PrimalStateVector&>(u).getVector();
1021  ROL::Ptr<const std::vector<Real> > zp =
1022  dynamic_cast<const PrimalControlVector&>(z).getVector();
1023  fem_->apply_adjoint_pde_control_hessian(*ahwvp,*wp,*vp,*up,*zp);
1024  }
1025  else {
1026  ahwv.zero();
1027  }
1028  }
1029  void applyAdjointHessian_22(ROL::Vector<Real> &ahwv, const ROL::Vector<Real> &w, const ROL::Vector<Real> &v,
1030  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
1031  if ( useHessian_ ) {
1032  ROL::Ptr<std::vector<Real> > ahwvp =
1033  dynamic_cast<DualControlVector&>(ahwv).getVector();
1034  ROL::Ptr<const std::vector<Real> > wp =
1035  dynamic_cast<const DualConstraintVector&>(w).getVector();
1036  ROL::Ptr<const std::vector<Real> > vp =
1037  dynamic_cast<const PrimalControlVector&>(v).getVector();
1038  ROL::Ptr<const std::vector<Real> > up =
1039  dynamic_cast<const PrimalStateVector&>(u).getVector();
1040  ROL::Ptr<const std::vector<Real> > zp =
1041  dynamic_cast<const PrimalControlVector&>(z).getVector();
1042  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 
1054  typedef H1VectorPrimal<Real> PrimalStateVector;
1055  typedef H1VectorDual<Real> DualStateVector;
1056 
1057  typedef L2VectorPrimal<Real> PrimalControlVector;
1058  typedef L2VectorDual<Real> DualControlVector;
1059 
1060  Real alpha_; // Penalty Parameter
1061  ROL::Ptr<BurgersFEM<Real> > fem_;
1062  ROL::Ptr<ROL::Vector<Real> > ud_;
1063  ROL::Ptr<ROL::Vector<Real> > diff_;
1064 
1065 public:
1066  Objective_BurgersControl(const ROL::Ptr<BurgersFEM<Real> > &fem,
1067  const ROL::Ptr<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  ROL::Ptr<const std::vector<Real> > up =
1074  dynamic_cast<const PrimalStateVector&>(u).getVector();
1075  ROL::Ptr<const std::vector<Real> > zp =
1076  dynamic_cast<const PrimalControlVector&>(z).getVector();
1077  ROL::Ptr<const std::vector<Real> > udp =
1078  dynamic_cast<const L2VectorPrimal<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  ROL::Ptr<std::vector<Real> > gp =
1089  dynamic_cast<DualStateVector&>(g).getVector();
1090  ROL::Ptr<const std::vector<Real> > up =
1091  dynamic_cast<const PrimalStateVector&>(u).getVector();
1092  ROL::Ptr<const std::vector<Real> > udp =
1093  dynamic_cast<const L2VectorPrimal<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  ROL::Ptr<std::vector<Real> > gp =
1104  dynamic_cast<DualControlVector&>(g).getVector();
1105  ROL::Ptr<const std::vector<Real> > zp =
1106  dynamic_cast<const PrimalControlVector&>(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 
1114  void hessVec_11( ROL::Vector<Real> &hv, const ROL::Vector<Real> &v,
1115  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
1116  ROL::Ptr<std::vector<Real> > hvp =
1117  dynamic_cast<DualStateVector&>(hv).getVector();
1118  ROL::Ptr<const std::vector<Real> > vp =
1119  dynamic_cast<const PrimalStateVector&>(v).getVector();
1120 
1121  fem_->apply_mass(*hvp,*vp);
1122  }
1123 
1124  void hessVec_12( ROL::Vector<Real> &hv, const ROL::Vector<Real> &v,
1125  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
1126  hv.zero();
1127  }
1128 
1129  void hessVec_21( ROL::Vector<Real> &hv, const ROL::Vector<Real> &v,
1130  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
1131  hv.zero();
1132  }
1133 
1134  void hessVec_22( ROL::Vector<Real> &hv, const ROL::Vector<Real> &v,
1135  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
1136  ROL::Ptr<std::vector<Real> > hvp =
1137  dynamic_cast<DualControlVector&>(hv).getVector();
1138  ROL::Ptr<const std::vector<Real> > vp =
1139  dynamic_cast<const PrimalControlVector&>(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>
1149 class L2BoundConstraint : public ROL::BoundConstraint<Real> {
1150 private:
1151  int dim_;
1152  std::vector<Real> x_lo_;
1153  std::vector<Real> x_up_;
1154  Real min_diff_;
1155  Real scale_;
1156  ROL::Ptr<BurgersFEM<Real> > fem_;
1157  ROL::Ptr<ROL::Vector<Real> > l_;
1158  ROL::Ptr<ROL::Vector<Real> > u_;
1159 
1160  void cast_vector(ROL::Ptr<std::vector<Real> > &xvec,
1161  ROL::Vector<Real> &x) const {
1162  try {
1163  xvec = dynamic_cast<L2VectorPrimal<Real>&>(x).getVector();
1164  }
1165  catch (std::exception &e) {
1166  xvec = dynamic_cast<L2VectorDual<Real>&>(x).getVector();
1167  }
1168  }
1169 
1170  void cast_const_vector(ROL::Ptr<const std::vector<Real> > &xvec,
1171  const ROL::Vector<Real> &x) const {
1172  try {
1173  xvec = dynamic_cast<const L2VectorPrimal<Real>&>(x).getVector();
1174  }
1175  catch (std::exception &e) {
1176  xvec = dynamic_cast<const L2VectorDual<Real>&>(x).getVector();
1177  }
1178  }
1179 
1180  void axpy(std::vector<Real> &out, const Real a,
1181  const std::vector<Real> &x, const std::vector<Real> &y) const{
1182  out.resize(dim_,0.0);
1183  for (unsigned i = 0; i < dim_; i++) {
1184  out[i] = a*x[i] + y[i];
1185  }
1186  }
1187 
1188  void projection(std::vector<Real> &x) {
1189  for ( int i = 0; i < dim_; i++ ) {
1190  x[i] = std::max(x_lo_[i],std::min(x_up_[i],x[i]));
1191  }
1192  }
1193 
1194 public:
1195  L2BoundConstraint(std::vector<Real> &l, std::vector<Real> &u,
1196  const ROL::Ptr<BurgersFEM<Real> > &fem, Real scale = 1.0)
1197  : x_lo_(l), x_up_(u), scale_(scale), fem_(fem) {
1198  dim_ = x_lo_.size();
1199  for ( int i = 0; i < dim_; i++ ) {
1200  if ( i == 0 ) {
1201  min_diff_ = x_up_[i] - x_lo_[i];
1202  }
1203  else {
1204  min_diff_ = ( (min_diff_ < (x_up_[i] - x_lo_[i])) ? min_diff_ : (x_up_[i] - x_lo_[i]) );
1205  }
1206  }
1207  min_diff_ *= 0.5;
1208  l_ = ROL::makePtr<L2VectorPrimal<Real>>(
1209  ROL::makePtr<std::vector<Real>>(l), fem);
1210  u_ = ROL::makePtr<L2VectorPrimal<Real>>(
1211  ROL::makePtr<std::vector<Real>>(u), fem);
1212  }
1213 
1214  bool isFeasible( const ROL::Vector<Real> &x ) {
1215  ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1216  bool val = true;
1217  int cnt = 1;
1218  for ( int i = 0; i < dim_; i++ ) {
1219  if ( (*ex)[i] >= x_lo_[i] && (*ex)[i] <= x_up_[i] ) { cnt *= 1; }
1220  else { cnt *= 0; }
1221  }
1222  if ( cnt == 0 ) { val = false; }
1223  return val;
1224  }
1225 
1226  void project( ROL::Vector<Real> &x ) {
1227  ROL::Ptr<std::vector<Real> > ex; cast_vector(ex,x);
1228  projection(*ex);
1229  }
1230 
1231  void pruneLowerActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &x, Real eps) {
1232  ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1233  ROL::Ptr<std::vector<Real> > ev; cast_vector(ev,v);
1234  Real epsn = std::min(scale_*eps,min_diff_);
1235  for ( int i = 0; i < dim_; i++ ) {
1236  if ( ((*ex)[i] <= x_lo_[i]+epsn) ) {
1237  (*ev)[i] = 0.0;
1238  }
1239  }
1240  }
1241 
1242  void pruneUpperActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &x, Real eps) {
1243  ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1244  ROL::Ptr<std::vector<Real> > ev; cast_vector(ev,v);
1245  Real epsn = std::min(scale_*eps,min_diff_);
1246  for ( int i = 0; i < dim_; i++ ) {
1247  if ( ((*ex)[i] >= x_up_[i]-epsn) ) {
1248  (*ev)[i] = 0.0;
1249  }
1250  }
1251  }
1252 
1253  void pruneActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &x, Real eps) {
1254  ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1255  ROL::Ptr<std::vector<Real> > ev; cast_vector(ev,v);
1256  Real epsn = std::min(scale_*eps,min_diff_);
1257  for ( int i = 0; i < dim_; i++ ) {
1258  if ( ((*ex)[i] <= x_lo_[i]+epsn) ||
1259  ((*ex)[i] >= x_up_[i]-epsn) ) {
1260  (*ev)[i] = 0.0;
1261  }
1262  }
1263  }
1264 
1265  void pruneLowerActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &g, const ROL::Vector<Real> &x, Real eps) {
1266  ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1267  ROL::Ptr<const std::vector<Real> > eg; cast_const_vector(eg,g);
1268  ROL::Ptr<std::vector<Real> > ev; cast_vector(ev,v);
1269  Real epsn = std::min(scale_*eps,min_diff_);
1270  for ( int i = 0; i < dim_; i++ ) {
1271  if ( ((*ex)[i] <= x_lo_[i]+epsn && (*eg)[i] > 0.0) ) {
1272  (*ev)[i] = 0.0;
1273  }
1274  }
1275  }
1276 
1277  void pruneUpperActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &g, const ROL::Vector<Real> &x, Real eps) {
1278  ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1279  ROL::Ptr<const std::vector<Real> > eg; cast_const_vector(eg,g);
1280  ROL::Ptr<std::vector<Real> > ev; cast_vector(ev,v);
1281  Real epsn = std::min(scale_*eps,min_diff_);
1282  for ( int i = 0; i < dim_; i++ ) {
1283  if ( ((*ex)[i] >= x_up_[i]-epsn && (*eg)[i] < 0.0) ) {
1284  (*ev)[i] = 0.0;
1285  }
1286  }
1287  }
1288 
1289  void pruneActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &g, const ROL::Vector<Real> &x, Real eps) {
1290  ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1291  ROL::Ptr<const std::vector<Real> > eg; cast_const_vector(eg,g);
1292  ROL::Ptr<std::vector<Real> > ev; cast_vector(ev,v);
1293  Real epsn = std::min(scale_*eps,min_diff_);
1294  for ( int i = 0; i < dim_; i++ ) {
1295  if ( ((*ex)[i] <= x_lo_[i]+epsn && (*eg)[i] > 0.0) ||
1296  ((*ex)[i] >= x_up_[i]-epsn && (*eg)[i] < 0.0) ) {
1297  (*ev)[i] = 0.0;
1298  }
1299  }
1300  }
1301 
1302  const ROL::Ptr<const ROL::Vector<Real> > getLowerBound(void) const {
1303  return l_;
1304  }
1305 
1306  const ROL::Ptr<const ROL::Vector<Real> > getUpperBound(void) const {
1307  return u_;
1308  }
1309 };
1310 
1311 template<class Real>
1312 class H1BoundConstraint : public ROL::BoundConstraint<Real> {
1313 private:
1314  int dim_;
1315  std::vector<Real> x_lo_;
1316  std::vector<Real> x_up_;
1317  Real min_diff_;
1318  Real scale_;
1319  ROL::Ptr<BurgersFEM<Real> > fem_;
1320  ROL::Ptr<ROL::Vector<Real> > l_;
1321  ROL::Ptr<ROL::Vector<Real> > u_;
1322 
1323  void cast_vector(ROL::Ptr<std::vector<Real> > &xvec,
1324  ROL::Vector<Real> &x) const {
1325  try {
1326  xvec = dynamic_cast<H1VectorPrimal<Real>&>(x).getVector();
1327  }
1328  catch (std::exception &e) {
1329  xvec = dynamic_cast<H1VectorDual<Real>&>(x).getVector();
1330  }
1331  }
1332 
1333  void cast_const_vector(ROL::Ptr<const std::vector<Real> > &xvec,
1334  const ROL::Vector<Real> &x) const {
1335  try {
1336  xvec = dynamic_cast<const H1VectorPrimal<Real>&>(x).getVector();
1337  }
1338  catch (std::exception &e) {
1339  xvec = dynamic_cast<const H1VectorDual<Real>&>(x).getVector();
1340  }
1341  }
1342 
1343  void axpy(std::vector<Real> &out, const Real a,
1344  const std::vector<Real> &x, const std::vector<Real> &y) const{
1345  out.resize(dim_,0.0);
1346  for (unsigned i = 0; i < dim_; i++) {
1347  out[i] = a*x[i] + y[i];
1348  }
1349  }
1350 
1351  void projection(std::vector<Real> &x) {
1352  for ( int i = 0; i < dim_; i++ ) {
1353  x[i] = std::max(x_lo_[i],std::min(x_up_[i],x[i]));
1354  }
1355  }
1356 
1357 public:
1358  H1BoundConstraint(std::vector<Real> &l, std::vector<Real> &u,
1359  const ROL::Ptr<BurgersFEM<Real> > &fem, Real scale = 1.0)
1360  : x_lo_(l), x_up_(u), scale_(scale), fem_(fem) {
1361  dim_ = x_lo_.size();
1362  for ( int i = 0; i < dim_; i++ ) {
1363  if ( i == 0 ) {
1364  min_diff_ = x_up_[i] - x_lo_[i];
1365  }
1366  else {
1367  min_diff_ = ( (min_diff_ < (x_up_[i] - x_lo_[i])) ? min_diff_ : (x_up_[i] - x_lo_[i]) );
1368  }
1369  }
1370  min_diff_ *= 0.5;
1371  l_ = ROL::makePtr<H1VectorPrimal<Real>>(
1372  ROL::makePtr<std::vector<Real>>(l), fem);
1373  u_ = ROL::makePtr<H1VectorPrimal<Real>>(
1374  ROL::makePtr<std::vector<Real>>(u), fem);
1375  }
1376 
1377  bool isFeasible( const ROL::Vector<Real> &x ) {
1378  ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1379  bool val = true;
1380  int cnt = 1;
1381  for ( int i = 0; i < dim_; i++ ) {
1382  if ( (*ex)[i] >= x_lo_[i] && (*ex)[i] <= x_up_[i] ) { cnt *= 1; }
1383  else { cnt *= 0; }
1384  }
1385  if ( cnt == 0 ) { val = false; }
1386  return val;
1387  }
1388 
1389  void project( ROL::Vector<Real> &x ) {
1390  ROL::Ptr<std::vector<Real> > ex; cast_vector(ex,x);
1391  projection(*ex);
1392  }
1393 
1394  void pruneLowerActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &x, Real eps) {
1395  ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1396  ROL::Ptr<std::vector<Real> > ev; cast_vector(ev,v);
1397  Real epsn = std::min(scale_*eps,min_diff_);
1398  for ( int i = 0; i < dim_; i++ ) {
1399  if ( ((*ex)[i] <= x_lo_[i]+epsn) ) {
1400  (*ev)[i] = 0.0;
1401  }
1402  }
1403  }
1404 
1405  void pruneUpperActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &x, Real eps) {
1406  ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1407  ROL::Ptr<std::vector<Real> > ev; cast_vector(ev,v);
1408  Real epsn = std::min(scale_*eps,min_diff_);
1409  for ( int i = 0; i < dim_; i++ ) {
1410  if ( ((*ex)[i] >= x_up_[i]-epsn) ) {
1411  (*ev)[i] = 0.0;
1412  }
1413  }
1414  }
1415 
1416  void pruneActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &x, Real eps) {
1417  ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1418  ROL::Ptr<std::vector<Real> > ev; cast_vector(ev,v);
1419  Real epsn = std::min(scale_*eps,min_diff_);
1420  for ( int i = 0; i < dim_; i++ ) {
1421  if ( ((*ex)[i] <= x_lo_[i]+epsn) ||
1422  ((*ex)[i] >= x_up_[i]-epsn) ) {
1423  (*ev)[i] = 0.0;
1424  }
1425  }
1426  }
1427 
1428  void pruneLowerActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &g, const ROL::Vector<Real> &x, Real eps) {
1429  ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1430  ROL::Ptr<const std::vector<Real> > eg; cast_const_vector(eg,g);
1431  ROL::Ptr<std::vector<Real> > ev; cast_vector(ev,v);
1432  Real epsn = std::min(scale_*eps,min_diff_);
1433  for ( int i = 0; i < dim_; i++ ) {
1434  if ( ((*ex)[i] <= x_lo_[i]+epsn && (*eg)[i] > 0.0) ) {
1435  (*ev)[i] = 0.0;
1436  }
1437  }
1438  }
1439 
1440  void pruneUpperActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &g, const ROL::Vector<Real> &x, Real eps) {
1441  ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1442  ROL::Ptr<const std::vector<Real> > eg; cast_const_vector(eg,g);
1443  ROL::Ptr<std::vector<Real> > ev; cast_vector(ev,v);
1444  Real epsn = std::min(scale_*eps,min_diff_);
1445  for ( int i = 0; i < dim_; i++ ) {
1446  if ( ((*ex)[i] >= x_up_[i]-epsn && (*eg)[i] < 0.0) ) {
1447  (*ev)[i] = 0.0;
1448  }
1449  }
1450  }
1451 
1452  void pruneActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &g, const ROL::Vector<Real> &x, Real eps) {
1453  ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1454  ROL::Ptr<const std::vector<Real> > eg; cast_const_vector(eg,g);
1455  ROL::Ptr<std::vector<Real> > ev; cast_vector(ev,v);
1456  Real epsn = std::min(scale_*eps,min_diff_);
1457  for ( int i = 0; i < dim_; i++ ) {
1458  if ( ((*ex)[i] <= x_lo_[i]+epsn && (*eg)[i] > 0.0) ||
1459  ((*ex)[i] >= x_up_[i]-epsn && (*eg)[i] < 0.0) ) {
1460  (*ev)[i] = 0.0;
1461  }
1462  }
1463  }
1464 
1465  const ROL::Ptr<const ROL::Vector<Real> > getLowerBound(void) const {
1466  return l_;
1467  }
1468 
1469  const ROL::Ptr<const ROL::Vector<Real> > getUpperBound(void) const {
1470  return u_;
1471  }
1472 };
Constraint_BurgersControl::DualConstraintVector
H1VectorPrimal< Real > DualConstraintVector
Definition: example_04.hpp:871
L2VectorDual::L2VectorDual
L2VectorDual(const ROL::Ptr< std::vector< Real > > &vec, const ROL::Ptr< BurgersFEM< Real > > &fem)
Definition: example_04.hpp:615
ROL::Objective_SimOpt
Provides the interface to evaluate simulation-based objective functions.
Definition: ROL_Objective_SimOpt.hpp:59
L2BoundConstraint::pruneUpperActive
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.
Definition: example_04.hpp:1242
L2BoundConstraint::x_up_
std::vector< Real > x_up_
Definition: example_04.hpp:1153
BurgersFEM::dx_
Real dx_
Definition: test_04.hpp:72
H1VectorDual::basis
ROL::Ptr< ROL::Vector< Real > > basis(const int i) const
Return i-th basis vector.
Definition: example_04.hpp:839
Objective_BurgersControl::diff_
ROL::Ptr< ROL::Vector< Real > > diff_
Definition: example_04.hpp:1063
BurgersFEM::apply_pde_jacobian
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
L2VectorDual::basis
ROL::Ptr< ROL::Vector< Real > > basis(const int i) const
Return i-th basis vector.
Definition: example_04.hpp:672
Constraint_BurgersControl::applyAdjointJacobian_1
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
L2BoundConstraint::axpy
void axpy(std::vector< Real > &out, const Real a, const std::vector< Real > &x, const std::vector< Real > &y) const
Definition: example_04.hpp:1180
BurgersFEM::cL2_
Real cL2_
Definition: test_04.hpp:79
L2VectorDual::norm
Real norm() const
Returns where .
Definition: example_04.hpp:654
L2BoundConstraint::isFeasible
bool isFeasible(const ROL::Vector< Real > &x)
Check if the vector, v, is feasible.
Definition: example_04.hpp:1214
H1VectorPrimal::dot
Real dot(const ROL::Vector< Real > &x) const
Compute where .
Definition: example_04.hpp:728
BurgersFEM::compute_residual
void compute_residual(std::vector< Real > &r, const std::vector< Real > &u, const std::vector< Real > &z) const
Definition: example_04.hpp:314
L2BoundConstraint::getUpperBound
const ROL::Ptr< const ROL::Vector< Real > > getUpperBound(void) const
Return the ref count pointer to the upper bound vector.
Definition: example_04.hpp:1306
L2VectorDual::getVector
ROL::Ptr< const std::vector< Real > > getVector() const
Definition: test_04.hpp:671
L2VectorPrimal::basis
ROL::Ptr< ROL::Vector< Real > > basis(const int i) const
Return i-th basis vector.
Definition: example_04.hpp:585
L2VectorDual::dual
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
BurgersFEM::apply_adjoint_pde_control_hessian
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
L2VectorPrimal::getVector
ROL::Ptr< const std::vector< Real > > getVector() const
Definition: test_04.hpp:584
Constraint_BurgersControl::applyAdjointHessian_12
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...
Definition: example_04.hpp:991
H1VectorPrimal::plus
void plus(const ROL::Vector< Real > &x)
Compute , where .
Definition: example_04.hpp:712
BurgersFEM::axpy
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
H1BoundConstraint::x_up_
std::vector< Real > x_up_
Definition: example_04.hpp:1316
BurgersFEM::apply_inverse_H1
void apply_inverse_H1(std::vector< Real > &Mu, const std::vector< Real > &u) const
Definition: example_04.hpp:285
Constraint_BurgersControl
Definition: test_04.hpp:868
H1VectorPrimal::H1VectorPrimal
H1VectorPrimal(const ROL::Ptr< std::vector< Real > > &vec, const ROL::Ptr< BurgersFEM< Real > > &fem)
Definition: example_04.hpp:702
L2VectorDual::clone
ROL::Ptr< ROL::Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
Definition: example_04.hpp:660
BurgersFEM::u0_
Real u0_
Definition: test_04.hpp:75
H1BoundConstraint::pruneActive
void pruneActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &g, const ROL::Vector< Real > &x, Real eps)
Definition: example_04.hpp:1452
ROL_Types.hpp
Contains definitions of custom data types in ROL.
H1VectorDual::dimension
int dimension() const
Return dimension of the vector space.
Definition: example_04.hpp:846
Constraint_BurgersControl::applyAdjointHessian_11
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 ...
Definition: example_04.hpp:971
H1VectorDual::plus
void plus(const ROL::Vector< Real > &x)
Compute , where .
Definition: example_04.hpp:792
L2BoundConstraint::pruneLowerActive
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.
Definition: example_04.hpp:1231
L2BoundConstraint::L2BoundConstraint
L2BoundConstraint(std::vector< Real > &l, std::vector< Real > &u, const ROL::Ptr< BurgersFEM< Real > > &fem, Real scale=1.0)
Definition: example_04.hpp:1195
Constraint_BurgersControl::applyInverseAdjointJacobian_1
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
L2BoundConstraint::pruneUpperActive
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.
Definition: example_04.hpp:1277
H1VectorDual::H1VectorDual
H1VectorDual(const ROL::Ptr< std::vector< Real > > &vec, const ROL::Ptr< BurgersFEM< Real > > &fem)
Definition: example_04.hpp:782
L2VectorPrimal::dimension
int dimension() const
Return dimension of the vector space.
Definition: example_04.hpp:592
BurgersFEM::compute_pde_jacobian
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
H1BoundConstraint::project
void project(ROL::Vector< Real > &x)
Project optimization variables onto the bounds.
Definition: example_04.hpp:1389
BurgersFEM::compute_H1_dot
Real compute_H1_dot(const std::vector< Real > &x, const std::vector< Real > &y) const
Definition: example_04.hpp:241
BurgersFEM::apply_control_jacobian
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
ROL::Vector::zero
virtual void zero()
Set to zero vector.
Definition: ROL_Vector.hpp:167
BurgersFEM::scale
void scale(std::vector< Real > &u, const Real alpha=0.0) const
Definition: example_04.hpp:95
H1VectorDual::getVector
ROL::Ptr< std::vector< Real > > getVector()
Definition: example_04.hpp:835
ROL::Vector
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:80
Objective_BurgersControl::fem_
ROL::Ptr< BurgersFEM< Real > > fem_
Definition: example_04.hpp:1061
Objective_BurgersControl::DualStateVector
H1VectorDual< Real > DualStateVector
Definition: example_04.hpp:1055
Constraint_BurgersControl::PrimalConstraintVector
H1VectorDual< Real > PrimalConstraintVector
Definition: example_04.hpp:870
L2VectorPrimal::scale
void scale(const Real alpha)
Compute where .
Definition: example_04.hpp:554
H1VectorPrimal::getVector
ROL::Ptr< const std::vector< Real > > getVector() const
Definition: test_04.hpp:751
L2VectorPrimal::L2VectorPrimal
L2VectorPrimal(const ROL::Ptr< std::vector< Real > > &vec, const ROL::Ptr< BurgersFEM< Real > > &fem)
Definition: example_04.hpp:535
Constraint_BurgersControl::PrimalStateVector
H1VectorPrimal< Real > PrimalStateVector
Definition: example_04.hpp:864
Objective_BurgersControl::hessVec_22
void hessVec_22(ROL::Vector< Real > &hv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Definition: example_04.hpp:1134
L2VectorPrimal::clone
ROL::Ptr< ROL::Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
Definition: example_04.hpp:573
Constraint_BurgersControl::Constraint_BurgersControl
Constraint_BurgersControl(ROL::Ptr< BurgersFEM< Real > > &fem, bool useHessian=true)
Definition: example_04.hpp:877
Objective_BurgersControl::value
Real value(const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Compute value.
Definition: example_04.hpp:1072
H1BoundConstraint::pruneActive
void pruneActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &x, Real eps)
Definition: example_04.hpp:1416
L2VectorPrimal::plus
void plus(const ROL::Vector< Real > &x)
Compute , where .
Definition: example_04.hpp:545
BurgersFEM::apply_inverse_pde_jacobian
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
Objective_BurgersControl::gradient_1
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.
Definition: example_04.hpp:1087
H1VectorPrimal::dual
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
H1BoundConstraint::isFeasible
bool isFeasible(const ROL::Vector< Real > &x)
Check if the vector, v, is feasible.
Definition: example_04.hpp:1377
BurgersFEM::test_inverse_mass
void test_inverse_mass(std::ostream &outStream=std::cout)
Definition: example_04.hpp:204
BurgersFEM::apply_adjoint_pde_hessian
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
BurgersFEM::apply_mass
void apply_mass(std::vector< Real > &Mu, const std::vector< Real > &u) const
Definition: example_04.hpp:174
H1BoundConstraint::getLowerBound
const ROL::Ptr< const ROL::Vector< Real > > getLowerBound(void) const
Return the ref count pointer to the lower bound vector.
Definition: example_04.hpp:1465
BurgersFEM::apply_H1
void apply_H1(std::vector< Real > &Mu, const std::vector< Real > &u) const
Definition: example_04.hpp:266
BurgersFEM::nl_
Real nl_
Definition: test_04.hpp:74
Objective_BurgersControl::ud_
ROL::Ptr< ROL::Vector< Real > > ud_
Definition: example_04.hpp:1062
L2BoundConstraint::pruneActive
void pruneActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &x, Real eps)
Definition: example_04.hpp:1253
BurgersFEM::cH1_
Real cH1_
Definition: test_04.hpp:78
Constraint_BurgersControl::applyAdjointHessian_21
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...
Definition: example_04.hpp:1010
BurgersFEM::u1_
Real u1_
Definition: test_04.hpp:76
H1BoundConstraint::pruneLowerActive
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 -binding set.
Definition: example_04.hpp:1428
H1VectorPrimal::norm
Real norm() const
Returns where .
Definition: example_04.hpp:734
BurgersFEM::compute_L2_dot
Real compute_L2_dot(const std::vector< Real > &x, const std::vector< Real > &y) const
Definition: example_04.hpp:151
L2BoundConstraint::fem_
ROL::Ptr< BurgersFEM< Real > > fem_
Definition: example_04.hpp:1156
L2BoundConstraint::cast_vector
void cast_vector(ROL::Ptr< std::vector< Real > > &xvec, ROL::Vector< Real > &x) const
Definition: example_04.hpp:1160
L2VectorDual::scale
void scale(const Real alpha)
Compute where .
Definition: example_04.hpp:634
BurgersFEM::mesh_spacing
Real mesh_spacing(void) const
Definition: example_04.hpp:143
BurgersFEM::apply_adjoint_pde_jacobian
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
L2VectorDual::dimension
int dimension() const
Return dimension of the vector space.
Definition: example_04.hpp:679
L2VectorPrimal::dual
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
H1VectorPrimal::clone
ROL::Ptr< ROL::Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
Definition: example_04.hpp:740
H1VectorDual::norm
Real norm() const
Returns where .
Definition: example_04.hpp:821
L2VectorDual::dot
Real dot(const ROL::Vector< Real > &x) const
Compute where .
Definition: example_04.hpp:641
ROL_BoundConstraint.hpp
BurgersFEM::apply_adjoint_control_jacobian
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
L2VectorPrimal::norm
Real norm() const
Returns where .
Definition: example_04.hpp:567
Constraint_BurgersControl::DualStateVector
H1VectorDual< Real > DualStateVector
Definition: example_04.hpp:865
H1BoundConstraint::getUpperBound
const ROL::Ptr< const ROL::Vector< Real > > getUpperBound(void) const
Return the ref count pointer to the upper bound vector.
Definition: example_04.hpp:1469
Constraint_BurgersControl::applyJacobian_1
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
ROL_Constraint_SimOpt.hpp
H1BoundConstraint::pruneLowerActive
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.
Definition: example_04.hpp:1394
L2VectorPrimal::dot
Real dot(const ROL::Vector< Real > &x) const
Compute where .
Definition: example_04.hpp:561
Objective_BurgersControl::hessVec_21
void hessVec_21(ROL::Vector< Real > &hv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Definition: example_04.hpp:1129
BurgersFEM::compute_H1_norm
Real compute_H1_norm(const std::vector< Real > &r) const
Definition: example_04.hpp:261
L2BoundConstraint::x_lo_
std::vector< Real > x_lo_
Definition: example_04.hpp:1152
H1VectorPrimal::getVector
ROL::Ptr< std::vector< Real > > getVector()
Definition: example_04.hpp:748
Constraint_BurgersControl::PrimalControlVector
L2VectorPrimal< Real > PrimalControlVector
Definition: example_04.hpp:867
H1BoundConstraint::u_
ROL::Ptr< ROL::Vector< Real > > u_
Definition: example_04.hpp:1321
ROL_Vector.hpp
Objective_BurgersControl::hessVec_12
void hessVec_12(ROL::Vector< Real > &hv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Definition: example_04.hpp:1124
L2VectorPrimal::getVector
ROL::Ptr< std::vector< Real > > getVector()
Definition: example_04.hpp:581
Objective_BurgersControl
Definition: burgers-control/example_01.hpp:64
L2BoundConstraint::pruneLowerActive
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 -binding set.
Definition: example_04.hpp:1265
H1VectorPrimal::scale
void scale(const Real alpha)
Compute where .
Definition: example_04.hpp:721
H1BoundConstraint::cast_const_vector
void cast_const_vector(ROL::Ptr< const std::vector< Real > > &xvec, const ROL::Vector< Real > &x) const
Definition: example_04.hpp:1333
H1BoundConstraint::projection
void projection(std::vector< Real > &x)
Definition: example_04.hpp:1351
Objective_BurgersControl::gradient_2
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.
Definition: example_04.hpp:1102
L2BoundConstraint::l_
ROL::Ptr< ROL::Vector< Real > > l_
Definition: example_04.hpp:1157
BurgersFEM::apply_adjoint_control_hessian
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
Constraint_BurgersControl::applyInverseJacobian_1
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
H1VectorPrimal::basis
ROL::Ptr< ROL::Vector< Real > > basis(const int i) const
Return i-th basis vector.
Definition: example_04.hpp:752
ROL::BoundConstraint
Provides the interface to apply upper and lower bound constraints.
Definition: ROL_BoundConstraint.hpp:72
L2BoundConstraint::projection
void projection(std::vector< Real > &x)
Definition: example_04.hpp:1188
H1BoundConstraint::axpy
void axpy(std::vector< Real > &out, const Real a, const std::vector< Real > &x, const std::vector< Real > &y) const
Definition: example_04.hpp:1343
H1VectorDual::dual
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
BurgersFEM::nu_
Real nu_
Definition: test_04.hpp:73
H1BoundConstraint::l_
ROL::Ptr< ROL::Vector< Real > > l_
Definition: example_04.hpp:1320
BurgersFEM::apply_inverse_mass
void apply_inverse_mass(std::vector< Real > &Mu, const std::vector< Real > &u) const
Definition: example_04.hpp:191
Objective_BurgersControl::hessVec_11
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.
Definition: example_04.hpp:1114
Constraint_BurgersControl::applyJacobian_2
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
BurgersFEM::test_inverse_H1
void test_inverse_H1(std::ostream &outStream=std::cout)
Definition: example_04.hpp:293
H1VectorDual::clone
ROL::Ptr< ROL::Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
Definition: example_04.hpp:827
BurgersFEM::linear_solve
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
Objective_BurgersControl::PrimalControlVector
L2VectorPrimal< Real > PrimalControlVector
Definition: example_04.hpp:1057
H1BoundConstraint::pruneUpperActive
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.
Definition: example_04.hpp:1440
L2BoundConstraint::pruneActive
void pruneActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &g, const ROL::Vector< Real > &x, Real eps)
Definition: example_04.hpp:1289
Objective_BurgersControl::PrimalStateVector
H1VectorPrimal< Real > PrimalStateVector
Definition: example_04.hpp:1054
H1BoundConstraint::x_lo_
std::vector< Real > x_lo_
Definition: example_04.hpp:1315
BurgersFEM::apply_inverse_adjoint_pde_jacobian
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
Constraint_BurgersControl::DualControlVector
L2VectorDual< Real > DualControlVector
Definition: example_04.hpp:868
BurgersFEM::BurgersFEM
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
BurgersFEM::f_
Real f_
Definition: test_04.hpp:77
BurgersFEM::compute_L2_norm
Real compute_L2_norm(const std::vector< Real > &r) const
Definition: example_04.hpp:169
Constraint_BurgersControl::applyAdjointHessian_22
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 ...
Definition: example_04.hpp:1029
L2VectorDual::plus
void plus(const ROL::Vector< Real > &x)
Compute , where .
Definition: example_04.hpp:625
H1VectorPrimal::dimension
int dimension() const
Return dimension of the vector space.
Definition: example_04.hpp:759
L2BoundConstraint::project
void project(ROL::Vector< Real > &x)
Project optimization variables onto the bounds.
Definition: example_04.hpp:1226
H1BoundConstraint::fem_
ROL::Ptr< BurgersFEM< Real > > fem_
Definition: example_04.hpp:1319
L2BoundConstraint::u_
ROL::Ptr< ROL::Vector< Real > > u_
Definition: example_04.hpp:1158
L2BoundConstraint::getLowerBound
const ROL::Ptr< const ROL::Vector< Real > > getLowerBound(void) const
Return the ref count pointer to the lower bound vector.
Definition: example_04.hpp:1302
BurgersFEM::update
void update(std::vector< Real > &u, const std::vector< Real > &s, const Real alpha=1.0) const
Definition: example_04.hpp:83
BurgersFEM::apply_adjoint_control_pde_hessian
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
H1BoundConstraint::cast_vector
void cast_vector(ROL::Ptr< std::vector< Real > > &xvec, ROL::Vector< Real > &x) const
Definition: example_04.hpp:1323
Constraint_BurgersControl::applyAdjointJacobian_2
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
BurgersFEM::num_dof
int num_dof(void) const
Definition: example_04.hpp:139
ROL::Constraint_SimOpt
Defines the constraint operator interface for simulation-based optimization.
Definition: ROL_Constraint_SimOpt.hpp:100
Objective_BurgersControl::DualControlVector
L2VectorDual< Real > DualControlVector
Definition: example_04.hpp:1058
Constraint_BurgersControl::value
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
H1VectorDual::dot
Real dot(const ROL::Vector< Real > &x) const
Compute where .
Definition: example_04.hpp:808
Objective_BurgersControl::Objective_BurgersControl
Objective_BurgersControl(const ROL::Ptr< BurgersFEM< Real > > &fem, const ROL::Ptr< ROL::Vector< Real > > &ud, Real alpha=1.e-4)
Definition: example_04.hpp:1066
H1VectorDual::getVector
ROL::Ptr< const std::vector< Real > > getVector() const
Definition: test_04.hpp:838
H1BoundConstraint::pruneUpperActive
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.
Definition: example_04.hpp:1405
L2BoundConstraint::cast_const_vector
void cast_const_vector(ROL::Ptr< const std::vector< Real > > &xvec, const ROL::Vector< Real > &x) const
Definition: example_04.hpp:1170
L2VectorDual::getVector
ROL::Ptr< std::vector< Real > > getVector()
Definition: example_04.hpp:668
ROL_Objective_SimOpt.hpp
H1BoundConstraint::H1BoundConstraint
H1BoundConstraint(std::vector< Real > &l, std::vector< Real > &u, const ROL::Ptr< BurgersFEM< Real > > &fem, Real scale=1.0)
Definition: example_04.hpp:1358
H1VectorDual::scale
void scale(const Real alpha)
Compute where .
Definition: example_04.hpp:801
Generated by   doxygen 1.8.13