ROL
example_06.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Rapid Optimization Library (ROL) Package
4 //
5 // Copyright 2014 NTESS and the ROL contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
15 #include "ROL_Types.hpp"
16 #include "ROL_Vector.hpp"
17 #include "ROL_BoundConstraint.hpp"
19 #include "ROL_Objective_SimOpt.hpp"
20 #include "ROL_TeuchosBatchManager.hpp"
21 
22 #include "Teuchos_LAPACK.hpp"
23 
24 template<class Real>
25 class L2VectorPrimal;
26 
27 template<class Real>
28 class L2VectorDual;
29 
30 template<class Real>
31 class H1VectorPrimal;
32 
33 template<class Real>
34 class H1VectorDual;
35 
36 template<class Real>
37 class BurgersFEM {
38 private:
39  int nx_;
40  Real dx_;
41  Real nu_;
42  Real nl_;
43  Real u0_;
44  Real u1_;
45  Real f_;
46  Real cH1_;
47  Real cL2_;
48 
49 private:
50  void update(std::vector<Real> &u, const std::vector<Real> &s, const Real alpha=1.0) const {
51  for (unsigned i=0; i<u.size(); i++) {
52  u[i] += alpha*s[i];
53  }
54  }
55 
56  void axpy(std::vector<Real> &out, const Real a, const std::vector<Real> &x, const std::vector<Real> &y) const {
57  for (unsigned i=0; i < x.size(); i++) {
58  out[i] = a*x[i] + y[i];
59  }
60  }
61 
62  void scale(std::vector<Real> &u, const Real alpha=0.0) const {
63  for (unsigned i=0; i<u.size(); i++) {
64  u[i] *= alpha;
65  }
66  }
67 
68  void linear_solve(std::vector<Real> &u, std::vector<Real> &dl, std::vector<Real> &d, std::vector<Real> &du,
69  const std::vector<Real> &r, const bool transpose = false) const {
70  if ( r.size() == 1 ) {
71  u.resize(1,r[0]/d[0]);
72  }
73  else if ( r.size() == 2 ) {
74  u.resize(2,0.0);
75  Real det = d[0]*d[1] - dl[0]*du[0];
76  u[0] = (d[1]*r[0] - du[0]*r[1])/det;
77  u[1] = (d[0]*r[1] - dl[0]*r[0])/det;
78  }
79  else {
80  u.assign(r.begin(),r.end());
81  // Perform LDL factorization
82  Teuchos::LAPACK<int,Real> lp;
83  std::vector<Real> du2(r.size()-2,0.0);
84  std::vector<int> ipiv(r.size(),0);
85  int info;
86  int dim = r.size();
87  int ldb = r.size();
88  int nhrs = 1;
89  lp.GTTRF(dim,&dl[0],&d[0],&du[0],&du2[0],&ipiv[0],&info);
90  char trans = 'N';
91  if ( transpose ) {
92  trans = 'T';
93  }
94  lp.GTTRS(trans,dim,nhrs,&dl[0],&d[0],&du[0],&du2[0],&ipiv[0],&u[0],ldb,&info);
95  }
96  }
97 
98 public:
99  BurgersFEM(int nx = 128, Real nl = 1.0, Real cH1 = 1.0, Real cL2 = 1.0)
100  : nx_(nx), dx_(1.0/((Real)nx+1.0)), nl_(nl), cH1_(cH1), cL2_(cL2) {}
101 
102  void set_problem_data(const Real nu, const Real f, const Real u0, const Real u1) {
103  nu_ = std::pow(10.0,nu-2.0);
104  f_ = 0.01*f;
105  u0_ = 1.0+0.001*u0;
106  u1_ = 0.001*u1;
107  }
108 
109  int num_dof(void) const {
110  return nx_;
111  }
112 
113  Real mesh_spacing(void) const {
114  return dx_;
115  }
116 
117  /***************************************************************************/
118  /*********************** L2 INNER PRODUCT FUNCTIONS ************************/
119  /***************************************************************************/
120  // Compute L2 inner product
121  Real compute_L2_dot(const std::vector<Real> &x, const std::vector<Real> &y) const {
122  Real ip = 0.0;
123  Real c = (((int)x.size()==nx_) ? 4.0 : 2.0);
124  for (unsigned i=0; i<x.size(); i++) {
125  if ( i == 0 ) {
126  ip += dx_/6.0*(c*x[i] + x[i+1])*y[i];
127  }
128  else if ( i == x.size()-1 ) {
129  ip += dx_/6.0*(x[i-1] + c*x[i])*y[i];
130  }
131  else {
132  ip += dx_/6.0*(x[i-1] + 4.0*x[i] + x[i+1])*y[i];
133  }
134  }
135  return ip;
136  }
137 
138  // compute L2 norm
139  Real compute_L2_norm(const std::vector<Real> &r) const {
140  return std::sqrt(compute_L2_dot(r,r));
141  }
142 
143  // Apply L2 Reisz operator
144  void apply_mass(std::vector<Real> &Mu, const std::vector<Real> &u ) const {
145  Mu.resize(u.size(),0.0);
146  Real c = (((int)u.size()==nx_) ? 4.0 : 2.0);
147  for (unsigned i=0; i<u.size(); i++) {
148  if ( i == 0 ) {
149  Mu[i] = dx_/6.0*(c*u[i] + u[i+1]);
150  }
151  else if ( i == u.size()-1 ) {
152  Mu[i] = dx_/6.0*(u[i-1] + c*u[i]);
153  }
154  else {
155  Mu[i] = dx_/6.0*(u[i-1] + 4.0*u[i] + u[i+1]);
156  }
157  }
158  }
159 
160  // Apply L2 inverse Reisz operator
161  void apply_inverse_mass(std::vector<Real> &Mu, const std::vector<Real> &u) const {
162  unsigned nx = u.size();
163  // Build mass matrix
164  std::vector<Real> dl(nx-1,dx_/6.0);
165  std::vector<Real> d(nx,2.0*dx_/3.0);
166  std::vector<Real> du(nx-1,dx_/6.0);
167  if ( (int)nx != nx_ ) {
168  d[ 0] = dx_/3.0;
169  d[nx-1] = dx_/3.0;
170  }
171  linear_solve(Mu,dl,d,du,u);
172  }
173 
174  void test_inverse_mass(std::ostream &outStream = std::cout) {
175  // State Mass Matrix
176  std::vector<Real> u(nx_,0.0), Mu(nx_,0.0), iMMu(nx_,0.0), diff(nx_,0.0);
177  for (int i = 0; i < nx_; i++) {
178  u[i] = 2.0*(Real)rand()/(Real)RAND_MAX - 1.0;
179  }
180  apply_mass(Mu,u);
181  apply_inverse_mass(iMMu,Mu);
182  axpy(diff,-1.0,iMMu,u);
183  Real error = compute_L2_norm(diff);
184  Real normu = compute_L2_norm(u);
185  outStream << "Test Inverse State Mass Matrix\n";
186  outStream << " ||u - inv(M)Mu|| = " << error << "\n";
187  outStream << " ||u|| = " << normu << "\n";
188  outStream << " Relative Error = " << error/normu << "\n";
189  outStream << "\n";
190  // Control Mass Matrix
191  u.resize(nx_+2,0.0); Mu.resize(nx_+2,0.0); iMMu.resize(nx_+2,0.0); diff.resize(nx_+2,0.0);
192  for (int i = 0; i < nx_+2; i++) {
193  u[i] = 2.0*(Real)rand()/(Real)RAND_MAX - 1.0;
194  }
195  apply_mass(Mu,u);
196  apply_inverse_mass(iMMu,Mu);
197  axpy(diff,-1.0,iMMu,u);
198  error = compute_L2_norm(diff);
199  normu = compute_L2_norm(u);
200  outStream << "Test Inverse Control Mass Matrix\n";
201  outStream << " ||z - inv(M)Mz|| = " << error << "\n";
202  outStream << " ||z|| = " << normu << "\n";
203  outStream << " Relative Error = " << error/normu << "\n";
204  outStream << "\n";
205  }
206 
207  /***************************************************************************/
208  /*********************** H1 INNER PRODUCT FUNCTIONS ************************/
209  /***************************************************************************/
210  // Compute H1 inner product
211  Real compute_H1_dot(const std::vector<Real> &x, const std::vector<Real> &y) const {
212  Real ip = 0.0;
213  for (int i=0; i<nx_; i++) {
214  if ( i == 0 ) {
215  ip += cL2_*dx_/6.0*(4.0*x[i] + x[i+1])*y[i]; // Mass term
216  ip += cH1_*(2.0*x[i] - x[i+1])/dx_*y[i]; // Stiffness term
217  }
218  else if ( i == nx_-1 ) {
219  ip += cL2_*dx_/6.0*(x[i-1] + 4.0*x[i])*y[i]; // Mass term
220  ip += cH1_*(2.0*x[i] - x[i-1])/dx_*y[i]; // Stiffness term
221  }
222  else {
223  ip += cL2_*dx_/6.0*(x[i-1] + 4.0*x[i] + x[i+1])*y[i]; // Mass term
224  ip += cH1_*(2.0*x[i] - x[i-1] - x[i+1])/dx_*y[i]; // Stiffness term
225  }
226  }
227  return ip;
228  }
229 
230  // compute H1 norm
231  Real compute_H1_norm(const std::vector<Real> &r) const {
232  return std::sqrt(compute_H1_dot(r,r));
233  }
234 
235  // Apply H2 Reisz operator
236  void apply_H1(std::vector<Real> &Mu, const std::vector<Real> &u ) const {
237  Mu.resize(nx_,0.0);
238  for (int i=0; i<nx_; i++) {
239  if ( i == 0 ) {
240  Mu[i] = cL2_*dx_/6.0*(4.0*u[i] + u[i+1])
241  + cH1_*(2.0*u[i] - u[i+1])/dx_;
242  }
243  else if ( i == nx_-1 ) {
244  Mu[i] = cL2_*dx_/6.0*(u[i-1] + 4.0*u[i])
245  + cH1_*(2.0*u[i] - u[i-1])/dx_;
246  }
247  else {
248  Mu[i] = cL2_*dx_/6.0*(u[i-1] + 4.0*u[i] + u[i+1])
249  + cH1_*(2.0*u[i] - u[i-1] - u[i+1])/dx_;
250  }
251  }
252  }
253 
254  // Apply H1 inverse Reisz operator
255  void apply_inverse_H1(std::vector<Real> &Mu, const std::vector<Real> &u) const {
256  // Build mass matrix
257  std::vector<Real> dl(nx_-1,cL2_*dx_/6.0 - cH1_/dx_);
258  std::vector<Real> d(nx_,2.0*(cL2_*dx_/3.0 + cH1_/dx_));
259  std::vector<Real> du(nx_-1,cL2_*dx_/6.0 - cH1_/dx_);
260  linear_solve(Mu,dl,d,du,u);
261  }
262 
263  void test_inverse_H1(std::ostream &outStream = std::cout) {
264  std::vector<Real> u(nx_,0.0), Mu(nx_,0.0), iMMu(nx_,0.0), diff(nx_,0.0);
265  for (int i = 0; i < nx_; i++) {
266  u[i] = 2.0*(Real)rand()/(Real)RAND_MAX - 1.0;
267  }
268  apply_H1(Mu,u);
269  apply_inverse_H1(iMMu,Mu);
270  axpy(diff,-1.0,iMMu,u);
271  Real error = compute_H1_norm(diff);
272  Real normu = compute_H1_norm(u);
273  outStream << "Test Inverse State H1 Matrix\n";
274  outStream << " ||u - inv(M)Mu|| = " << error << "\n";
275  outStream << " ||u|| = " << normu << "\n";
276  outStream << " Relative Error = " << error/normu << "\n";
277  outStream << "\n";
278  }
279 
280  /***************************************************************************/
281  /*********************** PDE RESIDUAL AND SOLVE ****************************/
282  /***************************************************************************/
283  // Compute current PDE residual
284  void compute_residual(std::vector<Real> &r, const std::vector<Real> &u,
285  const std::vector<Real> &z) const {
286  r.clear();
287  r.resize(nx_,0.0);
288  for (int i=0; i<nx_; i++) {
289  // Contribution from stiffness term
290  if ( i==0 ) {
291  r[i] = nu_/dx_*(2.0*u[i]-u[i+1]);
292  }
293  else if (i==nx_-1) {
294  r[i] = nu_/dx_*(2.0*u[i]-u[i-1]);
295  }
296  else {
297  r[i] = nu_/dx_*(2.0*u[i]-u[i-1]-u[i+1]);
298  }
299  // Contribution from nonlinear term
300  if (i<nx_-1){
301  r[i] += nl_*u[i+1]*(u[i]+u[i+1])/6.0;
302  }
303  if (i>0) {
304  r[i] -= nl_*u[i-1]*(u[i-1]+u[i])/6.0;
305  }
306  // Contribution from control
307  r[i] -= dx_/6.0*(z[i]+4.0*z[i+1]+z[i+2]);
308  // Contribution from right-hand side
309  r[i] -= dx_*f_;
310  }
311  // Contribution from Dirichlet boundary terms
312  r[0] -= nl_*(u0_*u[ 0]/6.0 + u0_*u0_/6.0) + nu_*u0_/dx_;
313  r[nx_-1] += nl_*(u1_*u[nx_-1]/6.0 + u1_*u1_/6.0) - nu_*u1_/dx_;
314  }
315 
316  /***************************************************************************/
317  /*********************** PDE JACOBIAN FUNCTIONS ****************************/
318  /***************************************************************************/
319  // Build PDE Jacobian trigiagonal matrix
320  void compute_pde_jacobian(std::vector<Real> &dl, // Lower diagonal
321  std::vector<Real> &d, // Diagonal
322  std::vector<Real> &du, // Upper diagonal
323  const std::vector<Real> &u) const { // State variable
324  // Get Diagonal and Off-Diagonal Entries of linear PDE Jacobian
325  d.clear();
326  d.resize(nx_,nu_*2.0/dx_);
327  dl.clear();
328  dl.resize(nx_-1,-nu_/dx_);
329  du.clear();
330  du.resize(nx_-1,-nu_/dx_);
331  // Contribution from nonlinearity
332  for (int i=0; i<nx_; i++) {
333  if (i<nx_-1) {
334  dl[i] += nl_*(-2.0*u[i]-u[i+1])/6.0;
335  d[i] += nl_*u[i+1]/6.0;
336  }
337  if (i>0) {
338  d[i] -= nl_*u[i-1]/6.0;
339  du[i-1] += nl_*(u[i-1]+2.0*u[i])/6.0;
340  }
341  }
342  // Contribution from Dirichlet boundary conditions
343  d[ 0] -= nl_*u0_/6.0;
344  d[nx_-1] += nl_*u1_/6.0;
345  }
346 
347  // Apply PDE Jacobian to a vector
348  void apply_pde_jacobian(std::vector<Real> &jv,
349  const std::vector<Real> &v,
350  const std::vector<Real> &u,
351  const std::vector<Real> &z) const {
352  // Fill jv
353  for (int i = 0; i < nx_; i++) {
354  jv[i] = nu_/dx_*2.0*v[i];
355  if ( i > 0 ) {
356  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]);
357  }
358  if ( i < nx_-1 ) {
359  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]);
360  }
361  }
362  jv[ 0] -= nl_*u0_/6.0*v[0];
363  jv[nx_-1] += nl_*u1_/6.0*v[nx_-1];
364  }
365 
366  // Apply inverse PDE Jacobian to a vector
367  void apply_inverse_pde_jacobian(std::vector<Real> &ijv,
368  const std::vector<Real> &v,
369  const std::vector<Real> &u,
370  const std::vector<Real> &z) const {
371  // Get PDE Jacobian
372  std::vector<Real> d(nx_,0.0);
373  std::vector<Real> dl(nx_-1,0.0);
374  std::vector<Real> du(nx_-1,0.0);
375  compute_pde_jacobian(dl,d,du,u);
376  // Solve solve state sensitivity system at current time step
377  linear_solve(ijv,dl,d,du,v);
378  }
379 
380  // Apply adjoint PDE Jacobian to a vector
381  void apply_adjoint_pde_jacobian(std::vector<Real> &ajv,
382  const std::vector<Real> &v,
383  const std::vector<Real> &u,
384  const std::vector<Real> &z) const {
385  // Fill jvp
386  for (int i = 0; i < nx_; i++) {
387  ajv[i] = nu_/dx_*2.0*v[i];
388  if ( i > 0 ) {
389  ajv[i] += -nu_/dx_*v[i-1]-nl_*(u[i-1]/6.0*v[i]
390  -(u[i-1]+2.0*u[i])/6.0*v[i-1]);
391  }
392  if ( i < nx_-1 ) {
393  ajv[i] += -nu_/dx_*v[i+1]+nl_*(u[i+1]/6.0*v[i]
394  -(u[i+1]+2.0*u[i])/6.0*v[i+1]);
395  }
396  }
397  ajv[ 0] -= nl_*u0_/6.0*v[0];
398  ajv[nx_-1] += nl_*u1_/6.0*v[nx_-1];
399  }
400 
401  // Apply inverse adjoint PDE Jacobian to a vector
402  void apply_inverse_adjoint_pde_jacobian(std::vector<Real> &iajv,
403  const std::vector<Real> &v,
404  const std::vector<Real> &u,
405  const std::vector<Real> &z) const {
406  // Get PDE Jacobian
407  std::vector<Real> d(nx_,0.0);
408  std::vector<Real> du(nx_-1,0.0);
409  std::vector<Real> dl(nx_-1,0.0);
410  compute_pde_jacobian(dl,d,du,u);
411  // Solve solve adjoint system at current time step
412  linear_solve(iajv,dl,d,du,v,true);
413  }
414 
415  /***************************************************************************/
416  /*********************** CONTROL JACOBIAN FUNCTIONS ************************/
417  /***************************************************************************/
418  // Apply control Jacobian to a vector
419  void apply_control_jacobian(std::vector<Real> &jv,
420  const std::vector<Real> &v,
421  const std::vector<Real> &u,
422  const std::vector<Real> &z) const {
423  for (int i=0; i<nx_; i++) {
424  // Contribution from control
425  jv[i] = -dx_/6.0*(v[i]+4.0*v[i+1]+v[i+2]);
426  }
427  }
428 
429  // Apply adjoint control Jacobian to a vector
430  void apply_adjoint_control_jacobian(std::vector<Real> &jv,
431  const std::vector<Real> &v,
432  const std::vector<Real> &u,
433  const std::vector<Real> &z) const {
434  for (int i=0; i<nx_+2; i++) {
435  if ( i == 0 ) {
436  jv[i] = -dx_/6.0*v[i];
437  }
438  else if ( i == 1 ) {
439  jv[i] = -dx_/6.0*(4.0*v[i-1]+v[i]);
440  }
441  else if ( i == nx_ ) {
442  jv[i] = -dx_/6.0*(4.0*v[i-1]+v[i-2]);
443  }
444  else if ( i == nx_+1 ) {
445  jv[i] = -dx_/6.0*v[i-2];
446  }
447  else {
448  jv[i] = -dx_/6.0*(v[i-2]+4.0*v[i-1]+v[i]);
449  }
450  }
451  }
452 
453  /***************************************************************************/
454  /*********************** AJDOINT HESSIANS **********************************/
455  /***************************************************************************/
456  void apply_adjoint_pde_hessian(std::vector<Real> &ahwv,
457  const std::vector<Real> &w,
458  const std::vector<Real> &v,
459  const std::vector<Real> &u,
460  const std::vector<Real> &z) const {
461  for (int i=0; i<nx_; i++) {
462  // Contribution from nonlinear term
463  ahwv[i] = 0.0;
464  if (i<nx_-1){
465  ahwv[i] += (w[i]*v[i+1] - w[i+1]*(2.0*v[i]+v[i+1]))/6.0;
466  }
467  if (i>0) {
468  ahwv[i] += (w[i-1]*(v[i-1]+2.0*v[i]) - w[i]*v[i-1])/6.0;
469  }
470  }
471  //ahwv.assign(u.size(),0.0);
472  }
473  void apply_adjoint_pde_control_hessian(std::vector<Real> &ahwv,
474  const std::vector<Real> &w,
475  const std::vector<Real> &v,
476  const std::vector<Real> &u,
477  const std::vector<Real> &z) {
478  ahwv.assign(u.size(),0.0);
479  }
480  void apply_adjoint_control_pde_hessian(std::vector<Real> &ahwv,
481  const std::vector<Real> &w,
482  const std::vector<Real> &v,
483  const std::vector<Real> &u,
484  const std::vector<Real> &z) {
485  ahwv.assign(z.size(),0.0);
486  }
487  void apply_adjoint_control_hessian(std::vector<Real> &ahwv,
488  const std::vector<Real> &w,
489  const std::vector<Real> &v,
490  const std::vector<Real> &u,
491  const std::vector<Real> &z) {
492  ahwv.assign(z.size(),0.0);
493  }
494 };
495 
496 template<class Real>
497 class L2VectorPrimal : public ROL::Vector<Real> {
498 private:
499  ROL::Ptr<std::vector<Real> > vec_;
500  ROL::Ptr<BurgersFEM<Real> > fem_;
501 
502  mutable ROL::Ptr<L2VectorDual<Real> > dual_vec_;
503 
504 public:
505  L2VectorPrimal(const ROL::Ptr<std::vector<Real> > & vec,
506  const ROL::Ptr<BurgersFEM<Real> > &fem)
507  : vec_(vec), fem_(fem), dual_vec_(ROL::nullPtr) {}
508 
509  void set( const ROL::Vector<Real> &x ) {
510  const L2VectorPrimal &ex = dynamic_cast<const L2VectorPrimal&>(x);
511  const std::vector<Real>& xval = *ex.getVector();
512  std::copy(xval.begin(),xval.end(),vec_->begin());
513  }
514 
515  void plus( const ROL::Vector<Real> &x ) {
516  const L2VectorPrimal &ex = dynamic_cast<const L2VectorPrimal&>(x);
517  const std::vector<Real>& xval = *ex.getVector();
518  unsigned dimension = vec_->size();
519  for (unsigned i=0; i<dimension; i++) {
520  (*vec_)[i] += xval[i];
521  }
522  }
523 
524  void scale( const Real alpha ) {
525  unsigned dimension = vec_->size();
526  for (unsigned i=0; i<dimension; i++) {
527  (*vec_)[i] *= alpha;
528  }
529  }
530 
531  Real dot( const ROL::Vector<Real> &x ) const {
532  const L2VectorPrimal & ex = dynamic_cast<const L2VectorPrimal&>(x);
533  const std::vector<Real>& xval = *ex.getVector();
534  return fem_->compute_L2_dot(xval,*vec_);
535  }
536 
537  Real norm() const {
538  Real val = 0;
539  val = std::sqrt( dot(*this) );
540  return val;
541  }
542 
543  ROL::Ptr<ROL::Vector<Real> > clone() const {
544  return ROL::makePtr<L2VectorPrimal>( ROL::makePtr<std::vector<Real>>(vec_->size(),0.0),fem_);
545  }
546 
547  ROL::Ptr<const std::vector<Real> > getVector() const {
548  return vec_;
549  }
550 
551  ROL::Ptr<std::vector<Real> > getVector() {
552  return vec_;
553  }
554 
555  ROL::Ptr<ROL::Vector<Real> > basis( const int i ) const {
556  ROL::Ptr<L2VectorPrimal> e
557  = ROL::makePtr<L2VectorPrimal>( ROL::makePtr<std::vector<Real>>(vec_->size(),0.0),fem_);
558  (*e->getVector())[i] = 1.0;
559  return e;
560  }
561 
562  int dimension() const {
563  return vec_->size();
564  }
565 
566  const ROL::Vector<Real>& dual() const {
567  dual_vec_ = ROL::makePtr<L2VectorDual<Real>>(
568  ROL::makePtr<std::vector<Real>>(*vec_),fem_);
569 
570  fem_->apply_mass(*(ROL::constPtrCast<std::vector<Real> >(dual_vec_->getVector())),*vec_);
571  return *dual_vec_;
572  }
573 
574  Real apply(const ROL::Vector<Real> &x) const {
575  const L2VectorDual<Real> &ex = dynamic_cast<const L2VectorDual<Real>&>(x);
576  const std::vector<Real>& xval = *ex.getVector();
577  return std::inner_product(vec_->begin(), vec_->end(), xval.begin(), Real(0));
578  }
579 
580 };
581 
582 template<class Real>
583 class L2VectorDual : public ROL::Vector<Real> {
584 private:
585  ROL::Ptr<std::vector<Real> > vec_;
586  ROL::Ptr<BurgersFEM<Real> > fem_;
587 
588  mutable ROL::Ptr<L2VectorPrimal<Real> > dual_vec_;
589 
590 public:
591  L2VectorDual(const ROL::Ptr<std::vector<Real> > & vec,
592  const ROL::Ptr<BurgersFEM<Real> > &fem)
593  : vec_(vec), fem_(fem), dual_vec_(ROL::nullPtr) {}
594 
595  void set( const ROL::Vector<Real> &x ) {
596  const L2VectorDual &ex = dynamic_cast<const L2VectorDual&>(x);
597  const std::vector<Real>& xval = *ex.getVector();
598  std::copy(xval.begin(),xval.end(),vec_->begin());
599  }
600 
601  void plus( const ROL::Vector<Real> &x ) {
602  const L2VectorDual &ex = dynamic_cast<const L2VectorDual&>(x);
603  const std::vector<Real>& xval = *ex.getVector();
604  unsigned dimension = vec_->size();
605  for (unsigned i=0; i<dimension; i++) {
606  (*vec_)[i] += xval[i];
607  }
608  }
609 
610  void scale( const Real alpha ) {
611  unsigned dimension = vec_->size();
612  for (unsigned i=0; i<dimension; i++) {
613  (*vec_)[i] *= alpha;
614  }
615  }
616 
617  Real dot( const ROL::Vector<Real> &x ) const {
618  const L2VectorDual & ex = dynamic_cast<const L2VectorDual&>(x);
619  const std::vector<Real>& xval = *ex.getVector();
620  unsigned dimension = vec_->size();
621  std::vector<Real> Mx(dimension,0.0);
622  fem_->apply_inverse_mass(Mx,xval);
623  Real val = 0.0;
624  for (unsigned i = 0; i < dimension; i++) {
625  val += Mx[i]*(*vec_)[i];
626  }
627  return val;
628  }
629 
630  Real norm() const {
631  Real val = 0;
632  val = std::sqrt( dot(*this) );
633  return val;
634  }
635 
636  ROL::Ptr<ROL::Vector<Real> > clone() const {
637  return ROL::makePtr<L2VectorDual>( ROL::makePtr<std::vector<Real>>(vec_->size(),0.0),fem_);
638  }
639 
640  ROL::Ptr<const std::vector<Real> > getVector() const {
641  return vec_;
642  }
643 
644  ROL::Ptr<std::vector<Real> > getVector() {
645  return vec_;
646  }
647 
648  ROL::Ptr<ROL::Vector<Real> > basis( const int i ) const {
649  ROL::Ptr<L2VectorDual> e
650  = ROL::makePtr<L2VectorDual>( ROL::makePtr<std::vector<Real>>(vec_->size(),0.0),fem_);
651  (*e->getVector())[i] = 1.0;
652  return e;
653  }
654 
655  int dimension() const {
656  return vec_->size();
657  }
658 
659  const ROL::Vector<Real>& dual() const {
660  dual_vec_ = ROL::makePtr<L2VectorPrimal<Real>>(
661  ROL::makePtr<std::vector<Real>>(*vec_),fem_);
662 
663  fem_->apply_inverse_mass(*(ROL::constPtrCast<std::vector<Real> >(dual_vec_->getVector())),*vec_);
664  return *dual_vec_;
665  }
666 
667  Real apply(const ROL::Vector<Real> &x) const {
668  const L2VectorPrimal<Real> &ex = dynamic_cast<const L2VectorPrimal<Real>&>(x);
669  const std::vector<Real>& xval = *ex.getVector();
670  return std::inner_product(vec_->begin(), vec_->end(), xval.begin(), Real(0));
671  }
672 
673 };
674 
675 template<class Real>
676 class H1VectorPrimal : public ROL::Vector<Real> {
677 private:
678  ROL::Ptr<std::vector<Real> > vec_;
679  ROL::Ptr<BurgersFEM<Real> > fem_;
680 
681  mutable ROL::Ptr<H1VectorDual<Real> > dual_vec_;
682 
683 public:
684  H1VectorPrimal(const ROL::Ptr<std::vector<Real> > & vec,
685  const ROL::Ptr<BurgersFEM<Real> > &fem)
686  : vec_(vec), fem_(fem), dual_vec_(ROL::nullPtr) {}
687 
688  void set( const ROL::Vector<Real> &x ) {
689  const H1VectorPrimal &ex = dynamic_cast<const H1VectorPrimal&>(x);
690  const std::vector<Real>& xval = *ex.getVector();
691  std::copy(xval.begin(),xval.end(),vec_->begin());
692  }
693 
694  void plus( const ROL::Vector<Real> &x ) {
695  const H1VectorPrimal &ex = dynamic_cast<const H1VectorPrimal&>(x);
696  const std::vector<Real>& xval = *ex.getVector();
697  unsigned dimension = vec_->size();
698  for (unsigned i=0; i<dimension; i++) {
699  (*vec_)[i] += xval[i];
700  }
701  }
702 
703  void scale( const Real alpha ) {
704  unsigned dimension = vec_->size();
705  for (unsigned i=0; i<dimension; i++) {
706  (*vec_)[i] *= alpha;
707  }
708  }
709 
710  Real dot( const ROL::Vector<Real> &x ) const {
711  const H1VectorPrimal & ex = dynamic_cast<const H1VectorPrimal&>(x);
712  const std::vector<Real>& xval = *ex.getVector();
713  return fem_->compute_H1_dot(xval,*vec_);
714  }
715 
716  Real norm() const {
717  Real val = 0;
718  val = std::sqrt( dot(*this) );
719  return val;
720  }
721 
722  ROL::Ptr<ROL::Vector<Real> > clone() const {
723  return ROL::makePtr<H1VectorPrimal>( ROL::makePtr<std::vector<Real>>(vec_->size(),0.0),fem_);
724  }
725 
726  ROL::Ptr<const std::vector<Real> > getVector() const {
727  return vec_;
728  }
729 
730  ROL::Ptr<std::vector<Real> > getVector() {
731  return vec_;
732  }
733 
734  ROL::Ptr<ROL::Vector<Real> > basis( const int i ) const {
735  ROL::Ptr<H1VectorPrimal> e
736  = ROL::makePtr<H1VectorPrimal>( ROL::makePtr<std::vector<Real>>(vec_->size(),0.0),fem_);
737  (*e->getVector())[i] = 1.0;
738  return e;
739  }
740 
741  int dimension() const {
742  return vec_->size();
743  }
744 
745  const ROL::Vector<Real>& dual() const {
746  dual_vec_ = ROL::makePtr<H1VectorDual<Real>>(
747  ROL::makePtr<std::vector<Real>>(*vec_),fem_);
748 
749  fem_->apply_H1(*(ROL::constPtrCast<std::vector<Real> >(dual_vec_->getVector())),*vec_);
750  return *dual_vec_;
751  }
752 
753  Real apply(const ROL::Vector<Real> &x) const {
754  const H1VectorDual<Real> &ex = dynamic_cast<const H1VectorDual<Real>&>(x);
755  const std::vector<Real>& xval = *ex.getVector();
756  return std::inner_product(vec_->begin(), vec_->end(), xval.begin(), Real(0));
757  }
758 
759 };
760 
761 template<class Real>
762 class H1VectorDual : public ROL::Vector<Real> {
763 private:
764  ROL::Ptr<std::vector<Real> > vec_;
765  ROL::Ptr<BurgersFEM<Real> > fem_;
766 
767  mutable ROL::Ptr<H1VectorPrimal<Real> > dual_vec_;
768 
769 public:
770  H1VectorDual(const ROL::Ptr<std::vector<Real> > & vec,
771  const ROL::Ptr<BurgersFEM<Real> > &fem)
772  : vec_(vec), fem_(fem), dual_vec_(ROL::nullPtr) {}
773 
774  void set( const ROL::Vector<Real> &x ) {
775  const H1VectorDual &ex = dynamic_cast<const H1VectorDual&>(x);
776  const std::vector<Real>& xval = *ex.getVector();
777  std::copy(xval.begin(),xval.end(),vec_->begin());
778  }
779 
780  void plus( const ROL::Vector<Real> &x ) {
781  const H1VectorDual &ex = dynamic_cast<const H1VectorDual&>(x);
782  const std::vector<Real>& xval = *ex.getVector();
783  unsigned dimension = vec_->size();
784  for (unsigned i=0; i<dimension; i++) {
785  (*vec_)[i] += xval[i];
786  }
787  }
788 
789  void scale( const Real alpha ) {
790  unsigned dimension = vec_->size();
791  for (unsigned i=0; i<dimension; i++) {
792  (*vec_)[i] *= alpha;
793  }
794  }
795 
796  Real dot( const ROL::Vector<Real> &x ) const {
797  const H1VectorDual & ex = dynamic_cast<const H1VectorDual&>(x);
798  const std::vector<Real>& xval = *ex.getVector();
799  unsigned dimension = vec_->size();
800  std::vector<Real> Mx(dimension,0.0);
801  fem_->apply_inverse_H1(Mx,xval);
802  Real val = 0.0;
803  for (unsigned i = 0; i < dimension; i++) {
804  val += Mx[i]*(*vec_)[i];
805  }
806  return val;
807  }
808 
809  Real norm() const {
810  Real val = 0;
811  val = std::sqrt( dot(*this) );
812  return val;
813  }
814 
815  ROL::Ptr<ROL::Vector<Real> > clone() const {
816  return ROL::makePtr<H1VectorDual>( ROL::makePtr<std::vector<Real>>(vec_->size(),0.0),fem_);
817  }
818 
819  ROL::Ptr<const std::vector<Real> > getVector() const {
820  return vec_;
821  }
822 
823  ROL::Ptr<std::vector<Real> > getVector() {
824  return vec_;
825  }
826 
827  ROL::Ptr<ROL::Vector<Real> > basis( const int i ) const {
828  ROL::Ptr<H1VectorDual> e
829  = ROL::makePtr<H1VectorDual>( ROL::makePtr<std::vector<Real>>(vec_->size(),0.0),fem_);
830  (*e->getVector())[i] = 1.0;
831  return e;
832  }
833 
834  int dimension() const {
835  return vec_->size();
836  }
837 
838  const ROL::Vector<Real>& dual() const {
839  dual_vec_ = ROL::makePtr<H1VectorPrimal<Real>>(
840  ROL::makePtr<std::vector<Real>>(*vec_),fem_);
841 
842  fem_->apply_inverse_H1(*(ROL::constPtrCast<std::vector<Real> >(dual_vec_->getVector())),*vec_);
843  return *dual_vec_;
844  }
845 
846  Real apply(const ROL::Vector<Real> &x) const {
847  const H1VectorPrimal<Real> &ex = dynamic_cast<const H1VectorPrimal<Real>&>(x);
848  const std::vector<Real>& xval = *ex.getVector();
849  return std::inner_product(vec_->begin(), vec_->end(), xval.begin(), Real(0));
850  }
851 
852 };
853 
854 template<class Real>
856 private:
857 
860 
863 
866 
867  ROL::Ptr<BurgersFEM<Real> > fem_;
868  bool useHessian_;
869 
870 public:
871  Constraint_BurgersControl(ROL::Ptr<BurgersFEM<Real> > &fem, bool useHessian = true)
872  : fem_(fem), useHessian_(useHessian) {}
873 
875  const ROL::Vector<Real> &z, Real &tol) {
876  ROL::Ptr<std::vector<Real> > cp =
877  dynamic_cast<PrimalConstraintVector&>(c).getVector();
878  ROL::Ptr<const std::vector<Real> > up =
879  dynamic_cast<const PrimalStateVector&>(u).getVector();
880  ROL::Ptr<const std::vector<Real> > zp =
881  dynamic_cast<const PrimalControlVector&>(z).getVector();
882 
883  const std::vector<Real> param
885  fem_->set_problem_data(param[0],param[1],param[2],param[3]);
886 
887  fem_->compute_residual(*cp,*up,*zp);
888  }
889 
891  const ROL::Vector<Real> &z, Real &tol) {
892  ROL::Ptr<std::vector<Real> > jvp =
893  dynamic_cast<PrimalConstraintVector&>(jv).getVector();
894  ROL::Ptr<const std::vector<Real> > vp =
895  dynamic_cast<const PrimalStateVector&>(v).getVector();
896  ROL::Ptr<const std::vector<Real> > up =
897  dynamic_cast<const PrimalStateVector&>(u).getVector();
898  ROL::Ptr<const std::vector<Real> > zp =
899  dynamic_cast<const PrimalControlVector&>(z).getVector();
900 
901  const std::vector<Real> param
903  fem_->set_problem_data(param[0],param[1],param[2],param[3]);
904 
905  fem_->apply_pde_jacobian(*jvp,*vp,*up,*zp);
906  }
907 
909  const ROL::Vector<Real> &z, Real &tol) {
910  ROL::Ptr<std::vector<Real> > jvp =
911  dynamic_cast<PrimalConstraintVector&>(jv).getVector();
912  ROL::Ptr<const std::vector<Real> > vp =
913  dynamic_cast<const PrimalControlVector&>(v).getVector();
914  ROL::Ptr<const std::vector<Real> > up =
915  dynamic_cast<const PrimalStateVector&>(u).getVector();
916  ROL::Ptr<const std::vector<Real> > zp =
917  dynamic_cast<const PrimalControlVector&>(z).getVector();
918 
919  const std::vector<Real> param
921  fem_->set_problem_data(param[0],param[1],param[2],param[3]);
922 
923  fem_->apply_control_jacobian(*jvp,*vp,*up,*zp);
924  }
925 
927  const ROL::Vector<Real> &z, Real &tol) {
928  ROL::Ptr<std::vector<Real> > ijvp =
929  dynamic_cast<PrimalStateVector&>(ijv).getVector();
930  ROL::Ptr<const std::vector<Real> > vp =
931  dynamic_cast<const PrimalConstraintVector&>(v).getVector();
932  ROL::Ptr<const std::vector<Real> > up =
933  dynamic_cast<const PrimalStateVector&>(u).getVector();
934  ROL::Ptr<const std::vector<Real> > zp =
935  dynamic_cast<const PrimalControlVector&>(z).getVector();
936 
937  const std::vector<Real> param
939  fem_->set_problem_data(param[0],param[1],param[2],param[3]);
940 
941  fem_->apply_inverse_pde_jacobian(*ijvp,*vp,*up,*zp);
942  }
943 
945  const ROL::Vector<Real> &z, Real &tol) {
946  ROL::Ptr<std::vector<Real> > jvp =
947  dynamic_cast<DualStateVector&>(ajv).getVector();
948  ROL::Ptr<const std::vector<Real> > vp =
949  dynamic_cast<const DualConstraintVector&>(v).getVector();
950  ROL::Ptr<const std::vector<Real> > up =
951  dynamic_cast<const PrimalStateVector&>(u).getVector();
952  ROL::Ptr<const std::vector<Real> > zp =
953  dynamic_cast<const PrimalControlVector&>(z).getVector();
954 
955  const std::vector<Real> param
957  fem_->set_problem_data(param[0],param[1],param[2],param[3]);
958 
959  fem_->apply_adjoint_pde_jacobian(*jvp,*vp,*up,*zp);
960  }
961 
963  const ROL::Vector<Real> &z, Real &tol) {
964  ROL::Ptr<std::vector<Real> > jvp =
965  dynamic_cast<DualControlVector&>(jv).getVector();
966  ROL::Ptr<const std::vector<Real> > vp =
967  dynamic_cast<const DualConstraintVector&>(v).getVector();
968  ROL::Ptr<const std::vector<Real> > up =
969  dynamic_cast<const PrimalStateVector&>(u).getVector();
970  ROL::Ptr<const std::vector<Real> > zp =
971  dynamic_cast<const PrimalControlVector&>(z).getVector();
972 
973  const std::vector<Real> param
975  fem_->set_problem_data(param[0],param[1],param[2],param[3]);
976 
977  fem_->apply_adjoint_control_jacobian(*jvp,*vp,*up,*zp);
978  }
979 
981  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
982  ROL::Ptr<std::vector<Real> > iajvp =
983  dynamic_cast<DualConstraintVector&>(iajv).getVector();
984  ROL::Ptr<const std::vector<Real> > vp =
985  dynamic_cast<const DualStateVector&>(v).getVector();
986  ROL::Ptr<const std::vector<Real> > up =
987  dynamic_cast<const PrimalStateVector&>(u).getVector();
988  ROL::Ptr<const std::vector<Real> > zp =
989  dynamic_cast<const PrimalControlVector&>(z).getVector();
990 
991  const std::vector<Real> param
993  fem_->set_problem_data(param[0],param[1],param[2],param[3]);
994 
995  fem_->apply_inverse_adjoint_pde_jacobian(*iajvp,*vp,*up,*zp);
996  }
997 
999  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
1000  if ( useHessian_ ) {
1001  ROL::Ptr<std::vector<Real> > ahwvp =
1002  dynamic_cast<DualStateVector&>(ahwv).getVector();
1003  ROL::Ptr<const std::vector<Real> > wp =
1004  dynamic_cast<const DualConstraintVector&>(w).getVector();
1005  ROL::Ptr<const std::vector<Real> > vp =
1006  dynamic_cast<const PrimalStateVector&>(v).getVector();
1007  ROL::Ptr<const std::vector<Real> > up =
1008  dynamic_cast<const PrimalStateVector&>(u).getVector();
1009  ROL::Ptr<const std::vector<Real> > zp =
1010  dynamic_cast<const PrimalControlVector&>(z).getVector();
1011 
1012  const std::vector<Real> param
1014  fem_->set_problem_data(param[0],param[1],param[2],param[3]);
1015 
1016  fem_->apply_adjoint_pde_hessian(*ahwvp,*wp,*vp,*up,*zp);
1017  }
1018  else {
1019  ahwv.zero();
1020  }
1021  }
1022 
1024  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
1025  if ( useHessian_ ) {
1026  ROL::Ptr<std::vector<Real> > ahwvp =
1027  dynamic_cast<DualControlVector&>(ahwv).getVector();
1028  ROL::Ptr<const std::vector<Real> > wp =
1029  dynamic_cast<const DualConstraintVector&>(w).getVector();
1030  ROL::Ptr<const std::vector<Real> > vp =
1031  dynamic_cast<const PrimalStateVector&>(v).getVector();
1032  ROL::Ptr<const std::vector<Real> > up =
1033  dynamic_cast<const PrimalStateVector&>(u).getVector();
1034  ROL::Ptr<const std::vector<Real> > zp =
1035  dynamic_cast<const PrimalControlVector&>(z).getVector();
1036 
1037  const std::vector<Real> param
1039  fem_->set_problem_data(param[0],param[1],param[2],param[3]);
1040 
1041  fem_->apply_adjoint_control_pde_hessian(*ahwvp,*wp,*vp,*up,*zp);
1042  }
1043  else {
1044  ahwv.zero();
1045  }
1046  }
1048  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
1049  if ( useHessian_ ) {
1050  ROL::Ptr<std::vector<Real> > ahwvp =
1051  dynamic_cast<DualStateVector&>(ahwv).getVector();
1052  ROL::Ptr<const std::vector<Real> > wp =
1053  dynamic_cast<const DualConstraintVector&>(w).getVector();
1054  ROL::Ptr<const std::vector<Real> > vp =
1055  dynamic_cast<const PrimalControlVector&>(v).getVector();
1056  ROL::Ptr<const std::vector<Real> > up =
1057  dynamic_cast<const PrimalStateVector&>(u).getVector();
1058  ROL::Ptr<const std::vector<Real> > zp =
1059  dynamic_cast<const PrimalControlVector&>(z).getVector();
1060 
1061  const std::vector<Real> param
1063  fem_->set_problem_data(param[0],param[1],param[2],param[3]);
1064 
1065  fem_->apply_adjoint_pde_control_hessian(*ahwvp,*wp,*vp,*up,*zp);
1066  }
1067  else {
1068  ahwv.zero();
1069  }
1070  }
1072  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
1073  if ( useHessian_ ) {
1074  ROL::Ptr<std::vector<Real> > ahwvp =
1075  dynamic_cast<DualControlVector&>(ahwv).getVector();
1076  ROL::Ptr<const std::vector<Real> > wp =
1077  dynamic_cast<const DualConstraintVector&>(w).getVector();
1078  ROL::Ptr<const std::vector<Real> > vp =
1079  dynamic_cast<const PrimalControlVector&>(v).getVector();
1080  ROL::Ptr<const std::vector<Real> > up =
1081  dynamic_cast<const PrimalStateVector&>(u).getVector();
1082  ROL::Ptr<const std::vector<Real> > zp =
1083  dynamic_cast<const PrimalControlVector&>(z).getVector();
1084 
1085  const std::vector<Real> param
1087  fem_->set_problem_data(param[0],param[1],param[2],param[3]);
1088 
1089  fem_->apply_adjoint_control_hessian(*ahwvp,*wp,*vp,*up,*zp);
1090  }
1091  else {
1092  ahwv.zero();
1093  }
1094  }
1095 };
1096 
1097 template<class Real>
1098 class Objective_BurgersControl : public ROL::Objective_SimOpt<Real> {
1099 private:
1100 
1103 
1106 
1107  Real alpha_; // Penalty Parameter
1108  ROL::Ptr<BurgersFEM<Real> > fem_;
1109  ROL::Ptr<ROL::Vector<Real> > ud_;
1110  ROL::Ptr<ROL::Vector<Real> > diff_;
1111 
1112 public:
1114  const ROL::Ptr<ROL::Vector<Real> > &ud,
1115  Real alpha = 1.e-4) : alpha_(alpha), fem_(fem), ud_(ud) {
1116  diff_ = ud_->clone();
1117  }
1118 
1120 
1121  Real value( const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
1122  ROL::Ptr<const std::vector<Real> > up =
1123  dynamic_cast<const PrimalStateVector&>(u).getVector();
1124  ROL::Ptr<const std::vector<Real> > zp =
1125  dynamic_cast<const PrimalControlVector&>(z).getVector();
1126  ROL::Ptr<const std::vector<Real> > udp =
1127  dynamic_cast<const L2VectorPrimal<Real>&>(*ud_).getVector();
1128 
1129  std::vector<Real> diff(udp->size(),0.0);
1130  for (unsigned i = 0; i < udp->size(); i++) {
1131  diff[i] = (*up)[i] - (*udp)[i];
1132  }
1133  return 0.5*(fem_->compute_L2_dot(diff,diff) + alpha_*fem_->compute_L2_dot(*zp,*zp));
1134  }
1135 
1136  void gradient_1( ROL::Vector<Real> &g, const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
1137  ROL::Ptr<std::vector<Real> > gp =
1138  dynamic_cast<DualStateVector&>(g).getVector();
1139  ROL::Ptr<const std::vector<Real> > up =
1140  dynamic_cast<const PrimalStateVector&>(u).getVector();
1141  ROL::Ptr<const std::vector<Real> > udp =
1142  dynamic_cast<const L2VectorPrimal<Real>&>(*ud_).getVector();
1143 
1144  std::vector<Real> diff(udp->size(),0.0);
1145  for (unsigned i = 0; i < udp->size(); i++) {
1146  diff[i] = (*up)[i] - (*udp)[i];
1147  }
1148  fem_->apply_mass(*gp,diff);
1149  }
1150 
1151  void gradient_2( ROL::Vector<Real> &g, const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
1152  ROL::Ptr<std::vector<Real> > gp =
1153  dynamic_cast<DualControlVector&>(g).getVector();
1154  ROL::Ptr<const std::vector<Real> > zp =
1155  dynamic_cast<const PrimalControlVector&>(z).getVector();
1156 
1157  fem_->apply_mass(*gp,*zp);
1158  for (unsigned i = 0; i < zp->size(); i++) {
1159  (*gp)[i] *= alpha_;
1160  }
1161  }
1162 
1164  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
1165  ROL::Ptr<std::vector<Real> > hvp =
1166  dynamic_cast<DualStateVector&>(hv).getVector();
1167  ROL::Ptr<const std::vector<Real> > vp =
1168  dynamic_cast<const PrimalStateVector&>(v).getVector();
1169 
1170  fem_->apply_mass(*hvp,*vp);
1171  }
1172 
1174  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
1175  hv.zero();
1176  }
1177 
1179  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
1180  hv.zero();
1181  }
1182 
1184  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
1185  ROL::Ptr<std::vector<Real> > hvp =
1186  dynamic_cast<DualControlVector&>(hv).getVector();
1187  ROL::Ptr<const std::vector<Real> > vp =
1188  dynamic_cast<const PrimalControlVector&>(v).getVector();
1189 
1190  fem_->apply_mass(*hvp,*vp);
1191  for (unsigned i = 0; i < vp->size(); i++) {
1192  (*hvp)[i] *= alpha_;
1193  }
1194  }
1195 };
1196 
1197 template<class Real, class Ordinal>
1198 class L2VectorBatchManager : public ROL::TeuchosBatchManager<Real,Ordinal> {
1199 private:
1200  void cast_vector(ROL::Ptr<std::vector<Real> > &xvec,
1201  ROL::Vector<Real> &x) const {
1202  try {
1203  xvec = dynamic_cast<L2VectorPrimal<Real>&>(x).getVector();
1204  }
1205  catch (std::exception &e) {
1206  xvec = dynamic_cast<L2VectorDual<Real>&>(x).getVector();
1207  }
1208  }
1209 
1210 public:
1211  L2VectorBatchManager(const ROL::Ptr<const Teuchos::Comm<int> > &comm)
1212  : ROL::TeuchosBatchManager<Real,Ordinal>(comm) {}
1214  ROL::Ptr<std::vector<Real> > input_ptr;
1215  cast_vector(input_ptr,input);
1216  int dim_i = input_ptr->size();
1217  ROL::Ptr<std::vector<Real> > output_ptr;
1218  cast_vector(output_ptr,output);
1219  int dim_o = output_ptr->size();
1220  if ( dim_i != dim_o ) {
1221  std::cout << "L2VectorBatchManager: DIMENSION MISMATCH ON RANK "
1222  << ROL::TeuchosBatchManager<Real,Ordinal>::batchID() << "\n";
1223  }
1224  else {
1225  ROL::TeuchosBatchManager<Real,Ordinal>::sumAll(&(*input_ptr)[0],&(*output_ptr)[0],dim_i);
1226  }
1227  }
1228 };
1229 
1230 template<class Real, class Ordinal>
1231 class H1VectorBatchManager : public ROL::TeuchosBatchManager<Real,Ordinal> {
1232 private:
1233  void cast_vector(ROL::Ptr<std::vector<Real> > &xvec,
1234  ROL::Vector<Real> &x) const {
1235  try {
1236  xvec = dynamic_cast<H1VectorPrimal<Real>&>(x).getVector();
1237  }
1238  catch (std::exception &e) {
1239  xvec = dynamic_cast<H1VectorDual<Real>&>(x).getVector();
1240  }
1241  }
1242 
1243 public:
1244  H1VectorBatchManager(const ROL::Ptr<const Teuchos::Comm<int> > &comm)
1245  : ROL::TeuchosBatchManager<Real,Ordinal>(comm) {}
1247  ROL::Ptr<std::vector<Real> > input_ptr;
1248  cast_vector(input_ptr,input);
1249  int dim_i = input_ptr->size();
1250  ROL::Ptr<std::vector<Real> > output_ptr;
1251  cast_vector(output_ptr,output);
1252  int dim_o = output_ptr->size();
1253  if ( dim_i != dim_o ) {
1254  std::cout << "H1VectorBatchManager: DIMENSION MISMATCH ON RANK "
1255  << ROL::TeuchosBatchManager<Real,Ordinal>::batchID() << "\n";
1256  }
1257  else {
1258  ROL::TeuchosBatchManager<Real,Ordinal>::sumAll(&(*input_ptr)[0],&(*output_ptr)[0],dim_i);
1259  }
1260  }
1261 };
1262 
1263 template<class Real>
1264 Real random(const ROL::Ptr<const Teuchos::Comm<int> > &comm) {
1265  Real val = 0.0;
1266  if ( Teuchos::rank<int>(*comm)==0 ) {
1267  val = (Real)rand()/(Real)RAND_MAX;
1268  }
1269  Teuchos::broadcast<int,Real>(*comm,0,1,&val);
1270  return val;
1271 }
H1VectorPrimal< Real > DualConstraintVector
Definition: example_06.hpp:865
BurgersFEM(int nx=128, Real nl=1.0, Real cH1=1.0, Real cL2=1.0)
Definition: example_06.hpp:99
L2VectorDual(const ROL::Ptr< std::vector< Real > > &vec, const ROL::Ptr< BurgersFEM< Real > > &fem)
Definition: example_06.hpp:591
ROL::Ptr< std::vector< Real > > vec_
Definition: example_04.hpp:674
Provides the interface to evaluate simulation-based objective functions.
Real dx_
Definition: test_04.hpp:38
ROL::Ptr< ROL::Vector< Real > > diff_
Real compute_H1_dot(const std::vector< Real > &x, const std::vector< Real > &y) const
Definition: example_06.hpp:211
void applyAdjointJacobian_1(ROL::Vector< Real > &ajv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Apply the adjoint of the partial constraint Jacobian at , , to the vector . This is the primary inter...
Definition: example_06.hpp:944
Real norm() const
Returns where .
Definition: example_06.hpp:716
Real cL2_
Definition: test_04.hpp:45
void apply_mass(std::vector< Real > &Mu, const std::vector< Real > &u) const
Definition: example_06.hpp:144
void axpy(std::vector< Real > &out, const Real a, const std::vector< Real > &x, const std::vector< Real > &y) const
Definition: example_06.hpp:56
ROL::Ptr< ROL::Vector< Real > > basis(const int i) const
Return i-th basis vector.
Definition: example_06.hpp:827
Real dot(const ROL::Vector< Real > &x) const override
Compute where .
Definition: test_04.hpp:636
Real dot(const ROL::Vector< Real > &x) const
Compute where .
Definition: example_06.hpp:617
void apply_adjoint_pde_control_hessian(std::vector< Real > &ahwv, const std::vector< Real > &w, const std::vector< Real > &v, const std::vector< Real > &u, const std::vector< Real > &z)
Definition: example_06.hpp:473
ROL::Ptr< std::vector< Real > > vec_
Definition: example_04.hpp:495
int dimension() const
Return dimension of the vector space.
Definition: example_04.hpp:558
ROL::Ptr< L2VectorDual< Real > > dual_vec_
Definition: test_04.hpp:503
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...
void plus(const ROL::Vector< Real > &x)
Compute , where .
Definition: example_06.hpp:694
void apply_adjoint_pde_hessian(std::vector< Real > &ahwv, const std::vector< Real > &w, const std::vector< Real > &v, const std::vector< Real > &u, const std::vector< Real > &z) const
Definition: example_06.hpp:456
ROL::Ptr< BurgersFEM< Real > > fem_
Definition: test_04.hpp:626
H1VectorPrimal(const ROL::Ptr< std::vector< Real > > &vec, const ROL::Ptr< BurgersFEM< Real > > &fem)
Definition: example_06.hpp:684
Real u0_
Definition: test_04.hpp:41
Real norm() const
Returns where .
Definition: example_06.hpp:537
Contains definitions of custom data types in ROL.
Real dot(const ROL::Vector< Real > &x) const
Compute where .
Definition: example_06.hpp:796
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_06.hpp:998
ROL::Ptr< const std::vector< Real > > getVector() const
Definition: example_04.hpp:543
void plus(const ROL::Vector< Real > &x)
Compute , where .
Definition: example_06.hpp:780
ROL::Ptr< BurgersFEM< Real > > fem_
Definition: test_04.hpp:587
const std::vector< Real > getParameter(void) const
Real compute_H1_norm(const std::vector< Real > &r) const
Definition: example_06.hpp:231
Real norm() const
Returns where .
Definition: example_06.hpp:630
Real apply(const ROL::Vector< Real > &x) const
Apply to a dual vector. This is equivalent to the call .
Definition: example_06.hpp:667
void applyInverseAdjointJacobian_1(ROL::Vector< Real > &iajv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Apply the inverse of the adjoint of the partial constraint Jacobian at , , to the vector ...
Definition: example_06.hpp:980
ROL::Ptr< std::vector< Real > > vec_
Definition: example_04.hpp:581
H1VectorDual(const ROL::Ptr< std::vector< Real > > &vec, const ROL::Ptr< BurgersFEM< Real > > &fem)
Definition: example_06.hpp:770
virtual void zero()
Set to zero vector.
Definition: ROL_Vector.hpp:133
ROL::Ptr< std::vector< Real > > getVector()
Definition: example_06.hpp:823
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:46
void cast_vector(ROL::Ptr< std::vector< Real > > &xvec, ROL::Vector< Real > &x) const
ROL::Ptr< BurgersFEM< Real > > fem_
H1VectorDual< Real > DualStateVector
H1VectorDual< Real > PrimalConstraintVector
Definition: example_06.hpp:864
Real apply(const ROL::Vector< Real > &x) const
Apply to a dual vector. This is equivalent to the call .
Definition: example_06.hpp:846
void scale(const Real alpha)
Compute where .
Definition: example_06.hpp:524
L2VectorPrimal(const ROL::Ptr< std::vector< Real > > &vec, const ROL::Ptr< BurgersFEM< Real > > &fem)
Definition: example_06.hpp:505
Real dot(const ROL::Vector< Real > &x) const override
Compute where .
Definition: test_04.hpp:597
Real dot(const ROL::Vector< Real > &x) const
Compute where .
Definition: example_06.hpp:710
int num_dof(void) const
Definition: example_06.hpp:109
H1VectorPrimal< Real > PrimalStateVector
Definition: example_06.hpp:858
void hessVec_22(ROL::Vector< Real > &hv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Constraint_BurgersControl(ROL::Ptr< BurgersFEM< Real > > &fem, bool useHessian=true)
Definition: example_06.hpp:871
Real value(const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Compute value.
void plus(const ROL::Vector< Real > &x)
Compute , where .
Definition: example_06.hpp:515
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_06.hpp:113
ROL::Ptr< ROL::Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
Definition: example_06.hpp:636
int dimension() const
Return dimension of the vector space.
Definition: example_04.hpp:737
ROL::Ptr< std::vector< Real > > vec_
Definition: example_04.hpp:760
void test_inverse_mass(std::ostream &outStream=std::cout)
Definition: example_06.hpp:174
void apply_pde_jacobian(std::vector< Real > &jv, const std::vector< Real > &v, const std::vector< Real > &u, const std::vector< Real > &z) const
Definition: example_06.hpp:348
ROL::Ptr< BurgersFEM< Real > > fem_
Definition: test_04.hpp:682
Real nl_
Definition: test_04.hpp:40
ROL::Ptr< ROL::Vector< Real > > ud_
L2VectorBatchManager(const ROL::Ptr< const Teuchos::Comm< int > > &comm)
ROL::Ptr< ROL::Vector< Real > > basis(const int i) const
Return i-th basis vector.
Definition: example_06.hpp:555
Real cH1_
Definition: test_04.hpp:44
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...
Real u1_
Definition: test_04.hpp:42
ROL::Ptr< const std::vector< Real > > getVector() const
Definition: example_04.hpp:722
ROL::Ptr< ROL::Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
Definition: example_06.hpp:815
ROL::Ptr< const std::vector< Real > > getVector() const
Definition: example_04.hpp:636
void set(const ROL::Vector< Real > &x)
Set where .
Definition: example_06.hpp:688
void scale(const Real alpha)
Compute where .
Definition: example_06.hpp:610
void compute_pde_jacobian(std::vector< Real > &dl, std::vector< Real > &d, std::vector< Real > &du, const std::vector< Real > &u) const
Definition: example_06.hpp:320
ROL::Ptr< ROL::Vector< Real > > basis(const int i) const
Return i-th basis vector.
Definition: example_06.hpp:648
H1VectorDual< Real > DualStateVector
Definition: example_06.hpp:859
void applyJacobian_1(ROL::Vector< Real > &jv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Apply the partial constraint Jacobian at , , to the vector .
Definition: example_06.hpp:890
Real compute_L2_dot(const std::vector< Real > &x, const std::vector< Real > &y) const
Definition: example_06.hpp:121
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_06.hpp:68
ROL::Ptr< std::vector< Real > > getVector()
Definition: example_06.hpp:730
L2VectorPrimal< Real > PrimalControlVector
Definition: example_06.hpp:861
void apply_inverse_mass(std::vector< Real > &Mu, const std::vector< Real > &u) const
Definition: example_06.hpp:161
void hessVec_12(ROL::Vector< Real > &hv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
ROL::Ptr< std::vector< Real > > getVector()
Definition: example_06.hpp:551
ROL::Ptr< BurgersFEM< Real > > fem_
Definition: test_04.hpp:502
ROL::Ptr< const std::vector< Real > > getVector() const
Definition: example_04.hpp:815
ROL::Ptr< BurgersFEM< Real > > fem_
Definition: test_04.hpp:541
void update(std::vector< Real > &u, const std::vector< Real > &s, const Real alpha=1.0) const
Definition: example_06.hpp:50
void scale(const Real alpha)
Compute where .
Definition: example_06.hpp:703
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.
Real random(const ROL::Ptr< const Teuchos::Comm< int > > &comm)
Definition: example_05.cpp:15
void set(const ROL::Vector< Real > &x)
Set where .
Definition: example_06.hpp:774
Real apply(const ROL::Vector< Real > &x) const
Apply to a dual vector. This is equivalent to the call .
Definition: example_06.hpp:753
void apply_adjoint_control_hessian(std::vector< Real > &ahwv, const std::vector< Real > &w, const std::vector< Real > &v, const std::vector< Real > &u, const std::vector< Real > &z)
Definition: example_06.hpp:487
void applyInverseJacobian_1(ROL::Vector< Real > &ijv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Apply the inverse partial constraint Jacobian at , , to the vector .
Definition: example_06.hpp:926
void compute_residual(std::vector< Real > &r, const std::vector< Real > &u, const std::vector< Real > &z) const
Definition: example_06.hpp:284
Real nu_
Definition: test_04.hpp:39
ROL::Ptr< ROL::Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
Definition: example_06.hpp:543
H1VectorBatchManager(const ROL::Ptr< const Teuchos::Comm< int > > &comm)
void set_problem_data(const Real nu, const Real f, const Real u0, const Real u1)
Definition: example_06.hpp:102
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 applyJacobian_2(ROL::Vector< Real > &jv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Apply the partial constraint Jacobian at , , to the vector .
Definition: example_06.hpp:908
void apply_control_jacobian(std::vector< Real > &jv, const std::vector< Real > &v, const std::vector< Real > &u, const std::vector< Real > &z) const
Definition: example_06.hpp:419
Real dot(const ROL::Vector< Real > &x) const override
Compute where .
Definition: test_04.hpp:551
void test_inverse_H1(std::ostream &outStream=std::cout)
Definition: example_06.hpp:263
const ROL::Vector< Real > & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
Definition: example_06.hpp:745
void apply_adjoint_control_jacobian(std::vector< Real > &jv, const std::vector< Real > &v, const std::vector< Real > &u, const std::vector< Real > &z) const
Definition: example_06.hpp:430
Real dot(const ROL::Vector< Real > &x) const
Compute where .
Definition: example_06.hpp:531
void cast_vector(ROL::Ptr< std::vector< Real > > &xvec, ROL::Vector< Real > &x) const
L2VectorPrimal< Real > PrimalControlVector
H1VectorPrimal< Real > PrimalStateVector
void apply_inverse_pde_jacobian(std::vector< Real > &ijv, const std::vector< Real > &v, const std::vector< Real > &u, const std::vector< Real > &z) const
Definition: example_06.hpp:367
L2VectorDual< Real > DualControlVector
Definition: example_06.hpp:862
Real f_
Definition: test_04.hpp:43
ROL::Ptr< ROL::Vector< Real > > basis(const int i) const
Return i-th basis vector.
Definition: example_06.hpp:734
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 ...
void apply_H1(std::vector< Real > &Mu, const std::vector< Real > &u) const
Definition: example_06.hpp:236
void plus(const ROL::Vector< Real > &x)
Compute , where .
Definition: example_06.hpp:601
void sumAll(ROL::Vector< Real > &input, ROL::Vector< Real > &output)
void apply_inverse_H1(std::vector< Real > &Mu, const std::vector< Real > &u) const
Definition: example_06.hpp:255
void apply_adjoint_pde_jacobian(std::vector< Real > &ajv, const std::vector< Real > &v, const std::vector< Real > &u, const std::vector< Real > &z) const
Definition: example_06.hpp:381
void set(const ROL::Vector< Real > &x)
Set where .
Definition: example_06.hpp:595
Real apply(const ROL::Vector< Real > &x) const
Apply to a dual vector. This is equivalent to the call .
Definition: example_06.hpp:574
ROL::Ptr< ROL::Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
Definition: example_06.hpp:722
void apply_inverse_adjoint_pde_jacobian(std::vector< Real > &iajv, const std::vector< Real > &v, const std::vector< Real > &u, const std::vector< Real > &z) const
Definition: example_06.hpp:402
void apply_adjoint_control_pde_hessian(std::vector< Real > &ahwv, const std::vector< Real > &w, const std::vector< Real > &v, const std::vector< Real > &u, const std::vector< Real > &z)
Definition: example_06.hpp:480
void applyAdjointJacobian_2(ROL::Vector< Real > &jv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Apply the adjoint of the partial constraint Jacobian at , , to vector . This is the primary interface...
Definition: example_06.hpp:962
int dimension() const
Return dimension of the vector space.
Definition: example_04.hpp:830
Defines the constraint operator interface for simulation-based optimization.
L2VectorDual< Real > DualControlVector
ROL::Ptr< H1VectorPrimal< Real > > dual_vec_
Definition: test_04.hpp:627
const ROL::Vector< Real > & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
Definition: example_06.hpp:838
void value(ROL::Vector< Real > &c, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Evaluate the constraint operator at .
Definition: example_06.hpp:874
void sumAll(ROL::Vector< Real > &input, ROL::Vector< Real > &output)
void set(const ROL::Vector< Real > &x)
Set where .
Definition: example_06.hpp:509
Real norm() const
Returns where .
Definition: example_06.hpp:809
void scale(std::vector< Real > &u, const Real alpha=0.0) const
Definition: example_06.hpp:62
const ROL::Vector< Real > & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
Definition: example_06.hpp:659
constexpr auto dim
ROL::Ptr< H1VectorDual< Real > > dual_vec_
Definition: test_04.hpp:588
Objective_BurgersControl(const ROL::Ptr< BurgersFEM< Real > > &fem, const ROL::Ptr< ROL::Vector< Real > > &ud, Real alpha=1.e-4)
Real dot(const ROL::Vector< Real > &x) const override
Compute where .
Definition: test_04.hpp:512
ROL::Ptr< std::vector< Real > > getVector()
Definition: example_06.hpp:644
void scale(const Real alpha)
Compute where .
Definition: example_06.hpp:789
int dimension() const
Return dimension of the vector space.
Definition: example_04.hpp:651
const ROL::Vector< Real > & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
Definition: example_06.hpp:566
Real compute_L2_norm(const std::vector< Real > &r) const
Definition: example_06.hpp:139
ROL::Ptr< L2VectorPrimal< Real > > dual_vec_
Definition: test_04.hpp:542