ROL
test_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 
48 #include "ROL_Types.hpp"
49 #include "ROL_Vector.hpp"
50 #include "ROL_BoundConstraint.hpp"
51 #include "ROL_Constraint_SimOpt.hpp"
52 #include "ROL_Objective_SimOpt.hpp"
53 
54 #include "Teuchos_LAPACK.hpp"
55 
56 template<class Real>
57 class L2VectorPrimal;
58 
59 template<class Real>
60 class L2VectorDual;
61 
62 template<class Real>
63 class H1VectorPrimal;
64 
65 template<class Real>
66 class H1VectorDual;
67 
68 template<class Real>
69 class BurgersFEM {
70 private:
71  int nx_;
72  Real dx_;
73  Real nu_;
74  Real nl_;
75  Real u0_;
76  Real u1_;
77  Real f_;
78  Real cH1_;
79  Real cL2_;
80 
81 private:
82  void update(std::vector<Real> &u, const std::vector<Real> &s, const Real alpha=1.0) const {
83  for (unsigned i=0; i<u.size(); i++) {
84  u[i] += alpha*s[i];
85  }
86  }
87 
88  void axpy(std::vector<Real> &out, const Real a, const std::vector<Real> &x, const std::vector<Real> &y) const {
89  for (unsigned i=0; i < x.size(); i++) {
90  out[i] = a*x[i] + y[i];
91  }
92  }
93 
94  void scale(std::vector<Real> &u, const Real alpha=0.0) const {
95  for (unsigned i=0; i<u.size(); i++) {
96  u[i] *= alpha;
97  }
98  }
99 
100  void linear_solve(std::vector<Real> &u, std::vector<Real> &dl, std::vector<Real> &d, std::vector<Real> &du,
101  const std::vector<Real> &r, const bool transpose = false) const {
102  if ( r.size() == 1 ) {
103  u.resize(1,r[0]/d[0]);
104  }
105  else if ( r.size() == 2 ) {
106  u.resize(2,0.0);
107  Real det = d[0]*d[1] - dl[0]*du[0];
108  u[0] = (d[1]*r[0] - du[0]*r[1])/det;
109  u[1] = (d[0]*r[1] - dl[0]*r[0])/det;
110  }
111  else {
112  u.assign(r.begin(),r.end());
113  // Perform LDL factorization
114  Teuchos::LAPACK<int,Real> lp;
115  std::vector<Real> du2(r.size()-2,0.0);
116  std::vector<int> ipiv(r.size(),0);
117  int info;
118  int dim = r.size();
119  int ldb = r.size();
120  int nhrs = 1;
121  lp.GTTRF(dim,&dl[0],&d[0],&du[0],&du2[0],&ipiv[0],&info);
122  char trans = 'N';
123  if ( transpose ) {
124  trans = 'T';
125  }
126  lp.GTTRS(trans,dim,nhrs,&dl[0],&d[0],&du[0],&du2[0],&ipiv[0],&u[0],ldb,&info);
127  }
128  }
129 
130 public:
131  BurgersFEM(int nx = 128, Real nl = 1.0, Real cH1 = 1.0, Real cL2 = 1.0)
132  : nx_(nx), dx_(1.0/((Real)nx+1.0)), nl_(nl), cH1_(cH1), cL2_(cL2) {
133  nu_ = 1.e-2;
134  f_ = 0.0;
135  u0_ = 1.0;
136  u1_ = 0.0;
137  }
138 
139  void set_problem_data(const Real nu, const Real f, const Real u0, const Real u1) {
140  nu_ = std::pow(10.0,nu-2.0);
141  f_ = 0.01*f;
142  u0_ = 1.0+0.001*u0;
143  u1_ = 0.001*u1;
144  }
145 
146  int num_dof(void) const {
147  return nx_;
148  }
149 
150  Real mesh_spacing(void) const {
151  return dx_;
152  }
153 
154  /***************************************************************************/
155  /*********************** L2 INNER PRODUCT FUNCTIONS ************************/
156  /***************************************************************************/
157  // Compute L2 inner product
158  Real compute_L2_dot(const std::vector<Real> &x, const std::vector<Real> &y) const {
159  Real ip = 0.0;
160  Real c = (((int)x.size()==nx_) ? 4.0 : 2.0);
161  for (unsigned i=0; i<x.size(); i++) {
162  if ( i == 0 ) {
163  ip += dx_/6.0*(c*x[i] + x[i+1])*y[i];
164  }
165  else if ( i == x.size()-1 ) {
166  ip += dx_/6.0*(x[i-1] + c*x[i])*y[i];
167  }
168  else {
169  ip += dx_/6.0*(x[i-1] + 4.0*x[i] + x[i+1])*y[i];
170  }
171  }
172  return ip;
173  }
174 
175  // compute L2 norm
176  Real compute_L2_norm(const std::vector<Real> &r) const {
177  return std::sqrt(compute_L2_dot(r,r));
178  }
179 
180  // Apply L2 Reisz operator
181  void apply_mass(std::vector<Real> &Mu, const std::vector<Real> &u ) const {
182  Mu.resize(u.size(),0.0);
183  Real c = (((int)u.size()==nx_) ? 4.0 : 2.0);
184  for (unsigned i=0; i<u.size(); i++) {
185  if ( i == 0 ) {
186  Mu[i] = dx_/6.0*(c*u[i] + u[i+1]);
187  }
188  else if ( i == u.size()-1 ) {
189  Mu[i] = dx_/6.0*(u[i-1] + c*u[i]);
190  }
191  else {
192  Mu[i] = dx_/6.0*(u[i-1] + 4.0*u[i] + u[i+1]);
193  }
194  }
195  }
196 
197  // Apply L2 inverse Reisz operator
198  void apply_inverse_mass(std::vector<Real> &Mu, const std::vector<Real> &u) const {
199  unsigned nx = u.size();
200  // Build mass matrix
201  std::vector<Real> dl(nx-1,dx_/6.0);
202  std::vector<Real> d(nx,2.0*dx_/3.0);
203  std::vector<Real> du(nx-1,dx_/6.0);
204  if ( (int)nx != nx_ ) {
205  d[ 0] = dx_/3.0;
206  d[nx-1] = dx_/3.0;
207  }
208  linear_solve(Mu,dl,d,du,u);
209  }
210 
211  void test_inverse_mass(std::ostream &outStream = std::cout) {
212  // State Mass Matrix
213  std::vector<Real> u(nx_,0.0), Mu(nx_,0.0), iMMu(nx_,0.0), diff(nx_,0.0);
214  for (int i = 0; i < nx_; i++) {
215  u[i] = 2.0*(Real)rand()/(Real)RAND_MAX - 1.0;
216  }
217  apply_mass(Mu,u);
218  apply_inverse_mass(iMMu,Mu);
219  axpy(diff,-1.0,iMMu,u);
220  Real error = compute_L2_norm(diff);
221  Real normu = compute_L2_norm(u);
222  outStream << "\nTest Inverse State Mass Matrix\n";
223  outStream << " ||u - inv(M)Mu|| = " << error << "\n";
224  outStream << " ||u|| = " << normu << "\n";
225  outStream << " Relative Error = " << error/normu << "\n";
226  outStream << "\n";
227  // Control Mass Matrix
228  u.resize(nx_+2,0.0); Mu.resize(nx_+2,0.0); iMMu.resize(nx_+2,0.0); diff.resize(nx_+2,0.0);
229  for (int i = 0; i < nx_+2; i++) {
230  u[i] = 2.0*(Real)rand()/(Real)RAND_MAX - 1.0;
231  }
232  apply_mass(Mu,u);
233  apply_inverse_mass(iMMu,Mu);
234  axpy(diff,-1.0,iMMu,u);
235  error = compute_L2_norm(diff);
236  normu = compute_L2_norm(u);
237  outStream << "\nTest Inverse Control Mass Matrix\n";
238  outStream << " ||z - inv(M)Mz|| = " << error << "\n";
239  outStream << " ||z|| = " << normu << "\n";
240  outStream << " Relative Error = " << error/normu << "\n";
241  outStream << "\n";
242  }
243 
244  /***************************************************************************/
245  /*********************** H1 INNER PRODUCT FUNCTIONS ************************/
246  /***************************************************************************/
247  // Compute H1 inner product
248  Real compute_H1_dot(const std::vector<Real> &x, const std::vector<Real> &y) const {
249  Real ip = 0.0;
250  for (int i=0; i<nx_; i++) {
251  if ( i == 0 ) {
252  ip += cL2_*dx_/6.0*(4.0*x[i] + x[i+1])*y[i]; // Mass term
253  ip += cH1_*(2.0*x[i] - x[i+1])/dx_*y[i]; // Stiffness term
254  }
255  else if ( i == nx_-1 ) {
256  ip += cL2_*dx_/6.0*(x[i-1] + 4.0*x[i])*y[i]; // Mass term
257  ip += cH1_*(2.0*x[i] - x[i-1])/dx_*y[i]; // Stiffness term
258  }
259  else {
260  ip += cL2_*dx_/6.0*(x[i-1] + 4.0*x[i] + x[i+1])*y[i]; // Mass term
261  ip += cH1_*(2.0*x[i] - x[i-1] - x[i+1])/dx_*y[i]; // Stiffness term
262  }
263  }
264  return ip;
265  }
266 
267  // compute H1 norm
268  Real compute_H1_norm(const std::vector<Real> &r) const {
269  return std::sqrt(compute_H1_dot(r,r));
270  }
271 
272  // Apply H2 Reisz operator
273  void apply_H1(std::vector<Real> &Mu, const std::vector<Real> &u ) const {
274  Mu.resize(nx_,0.0);
275  for (int i=0; i<nx_; i++) {
276  if ( i == 0 ) {
277  Mu[i] = cL2_*dx_/6.0*(4.0*u[i] + u[i+1])
278  + cH1_*(2.0*u[i] - u[i+1])/dx_;
279  }
280  else if ( i == nx_-1 ) {
281  Mu[i] = cL2_*dx_/6.0*(u[i-1] + 4.0*u[i])
282  + cH1_*(2.0*u[i] - u[i-1])/dx_;
283  }
284  else {
285  Mu[i] = cL2_*dx_/6.0*(u[i-1] + 4.0*u[i] + u[i+1])
286  + cH1_*(2.0*u[i] - u[i-1] - u[i+1])/dx_;
287  }
288  }
289  }
290 
291  // Apply H1 inverse Reisz operator
292  void apply_inverse_H1(std::vector<Real> &Mu, const std::vector<Real> &u) const {
293  // Build mass matrix
294  std::vector<Real> dl(nx_-1,cL2_*dx_/6.0 - cH1_/dx_);
295  std::vector<Real> d(nx_,2.0*(cL2_*dx_/3.0 + cH1_/dx_));
296  std::vector<Real> du(nx_-1,cL2_*dx_/6.0 - cH1_/dx_);
297  linear_solve(Mu,dl,d,du,u);
298  }
299 
300  void test_inverse_H1(std::ostream &outStream = std::cout) {
301  std::vector<Real> u(nx_,0.0), Mu(nx_,0.0), iMMu(nx_,0.0), diff(nx_,0.0);
302  for (int i = 0; i < nx_; i++) {
303  u[i] = 2.0*(Real)rand()/(Real)RAND_MAX - 1.0;
304  }
305  apply_H1(Mu,u);
306  apply_inverse_H1(iMMu,Mu);
307  axpy(diff,-1.0,iMMu,u);
308  Real error = compute_H1_norm(diff);
309  Real normu = compute_H1_norm(u);
310  outStream << "\nTest Inverse State H1 Matrix\n";
311  outStream << " ||u - inv(M)Mu|| = " << error << "\n";
312  outStream << " ||u|| = " << normu << "\n";
313  outStream << " Relative Error = " << error/normu << "\n";
314  outStream << "\n";
315  }
316 
317  /***************************************************************************/
318  /*********************** PDE RESIDUAL AND SOLVE ****************************/
319  /***************************************************************************/
320  // Compute current PDE residual
321  void compute_residual(std::vector<Real> &r, const std::vector<Real> &u,
322  const std::vector<Real> &z) const {
323  r.clear();
324  r.resize(nx_,0.0);
325  for (int i=0; i<nx_; i++) {
326  // Contribution from stiffness term
327  if ( i==0 ) {
328  r[i] = nu_/dx_*(2.0*u[i]-u[i+1]);
329  }
330  else if (i==nx_-1) {
331  r[i] = nu_/dx_*(2.0*u[i]-u[i-1]);
332  }
333  else {
334  r[i] = nu_/dx_*(2.0*u[i]-u[i-1]-u[i+1]);
335  }
336  // Contribution from nonlinear term
337  if (i<nx_-1){
338  r[i] += nl_*u[i+1]*(u[i]+u[i+1])/6.0;
339  }
340  if (i>0) {
341  r[i] -= nl_*u[i-1]*(u[i-1]+u[i])/6.0;
342  }
343  // Contribution from control
344  r[i] -= dx_/6.0*(z[i]+4.0*z[i+1]+z[i+2]);
345  // Contribution from right-hand side
346  r[i] -= dx_*f_;
347  }
348  // Contribution from Dirichlet boundary terms
349  r[0] -= nl_*(u0_*u[ 0]/6.0 + u0_*u0_/6.0) + nu_*u0_/dx_;
350  r[nx_-1] += nl_*(u1_*u[nx_-1]/6.0 + u1_*u1_/6.0) - nu_*u1_/dx_;
351  }
352 
353  /***************************************************************************/
354  /*********************** PDE JACOBIAN FUNCTIONS ****************************/
355  /***************************************************************************/
356  // Build PDE Jacobian trigiagonal matrix
357  void compute_pde_jacobian(std::vector<Real> &dl, // Lower diagonal
358  std::vector<Real> &d, // Diagonal
359  std::vector<Real> &du, // Upper diagonal
360  const std::vector<Real> &u) const { // State variable
361  // Get Diagonal and Off-Diagonal Entries of linear PDE Jacobian
362  d.clear();
363  d.resize(nx_,nu_*2.0/dx_);
364  dl.clear();
365  dl.resize(nx_-1,-nu_/dx_);
366  du.clear();
367  du.resize(nx_-1,-nu_/dx_);
368  // Contribution from nonlinearity
369  for (int i=0; i<nx_; i++) {
370  if (i<nx_-1) {
371  dl[i] += nl_*(-2.0*u[i]-u[i+1])/6.0;
372  d[i] += nl_*u[i+1]/6.0;
373  }
374  if (i>0) {
375  d[i] -= nl_*u[i-1]/6.0;
376  du[i-1] += nl_*(u[i-1]+2.0*u[i])/6.0;
377  }
378  }
379  // Contribution from Dirichlet boundary conditions
380  d[ 0] -= nl_*u0_/6.0;
381  d[nx_-1] += nl_*u1_/6.0;
382  }
383 
384  // Apply PDE Jacobian to a vector
385  void apply_pde_jacobian(std::vector<Real> &jv,
386  const std::vector<Real> &v,
387  const std::vector<Real> &u,
388  const std::vector<Real> &z) const {
389  // Fill jv
390  for (int i = 0; i < nx_; i++) {
391  jv[i] = nu_/dx_*2.0*v[i];
392  if ( i > 0 ) {
393  jv[i] += -nu_/dx_*v[i-1]-nl_*(u[i-1]/6.0*v[i]+(u[i]+2.0*u[i-1])/6.0*v[i-1]);
394  }
395  if ( i < nx_-1 ) {
396  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]);
397  }
398  }
399  jv[ 0] -= nl_*u0_/6.0*v[0];
400  jv[nx_-1] += nl_*u1_/6.0*v[nx_-1];
401  }
402 
403  // Apply inverse PDE Jacobian to a vector
404  void apply_inverse_pde_jacobian(std::vector<Real> &ijv,
405  const std::vector<Real> &v,
406  const std::vector<Real> &u,
407  const std::vector<Real> &z) const {
408  // Get PDE Jacobian
409  std::vector<Real> d(nx_,0.0);
410  std::vector<Real> dl(nx_-1,0.0);
411  std::vector<Real> du(nx_-1,0.0);
412  compute_pde_jacobian(dl,d,du,u);
413  // Solve solve state sensitivity system at current time step
414  linear_solve(ijv,dl,d,du,v);
415  }
416 
417  // Apply adjoint PDE Jacobian to a vector
418  void apply_adjoint_pde_jacobian(std::vector<Real> &ajv,
419  const std::vector<Real> &v,
420  const std::vector<Real> &u,
421  const std::vector<Real> &z) const {
422  // Fill jvp
423  for (int i = 0; i < nx_; i++) {
424  ajv[i] = nu_/dx_*2.0*v[i];
425  if ( i > 0 ) {
426  ajv[i] += -nu_/dx_*v[i-1]-nl_*(u[i-1]/6.0*v[i]
427  -(u[i-1]+2.0*u[i])/6.0*v[i-1]);
428  }
429  if ( i < nx_-1 ) {
430  ajv[i] += -nu_/dx_*v[i+1]+nl_*(u[i+1]/6.0*v[i]
431  -(u[i+1]+2.0*u[i])/6.0*v[i+1]);
432  }
433  }
434  ajv[ 0] -= nl_*u0_/6.0*v[0];
435  ajv[nx_-1] += nl_*u1_/6.0*v[nx_-1];
436  }
437 
438  // Apply inverse adjoint PDE Jacobian to a vector
439  void apply_inverse_adjoint_pde_jacobian(std::vector<Real> &iajv,
440  const std::vector<Real> &v,
441  const std::vector<Real> &u,
442  const std::vector<Real> &z) const {
443  // Get PDE Jacobian
444  std::vector<Real> d(nx_,0.0);
445  std::vector<Real> du(nx_-1,0.0);
446  std::vector<Real> dl(nx_-1,0.0);
447  compute_pde_jacobian(dl,d,du,u);
448  // Solve solve adjoint system at current time step
449  linear_solve(iajv,dl,d,du,v,true);
450  }
451 
452  /***************************************************************************/
453  /*********************** CONTROL JACOBIAN FUNCTIONS ************************/
454  /***************************************************************************/
455  // Apply control Jacobian to a vector
456  void apply_control_jacobian(std::vector<Real> &jv,
457  const std::vector<Real> &v,
458  const std::vector<Real> &u,
459  const std::vector<Real> &z) const {
460  for (int i=0; i<nx_; i++) {
461  // Contribution from control
462  jv[i] = -dx_/6.0*(v[i]+4.0*v[i+1]+v[i+2]);
463  }
464  }
465 
466  // Apply adjoint control Jacobian to a vector
467  void apply_adjoint_control_jacobian(std::vector<Real> &jv,
468  const std::vector<Real> &v,
469  const std::vector<Real> &u,
470  const std::vector<Real> &z) const {
471  for (int i=0; i<nx_+2; i++) {
472  if ( i == 0 ) {
473  jv[i] = -dx_/6.0*v[i];
474  }
475  else if ( i == 1 ) {
476  jv[i] = -dx_/6.0*(4.0*v[i-1]+v[i]);
477  }
478  else if ( i == nx_ ) {
479  jv[i] = -dx_/6.0*(4.0*v[i-1]+v[i-2]);
480  }
481  else if ( i == nx_+1 ) {
482  jv[i] = -dx_/6.0*v[i-2];
483  }
484  else {
485  jv[i] = -dx_/6.0*(v[i-2]+4.0*v[i-1]+v[i]);
486  }
487  }
488  }
489 
490  /***************************************************************************/
491  /*********************** AJDOINT HESSIANS **********************************/
492  /***************************************************************************/
493  void apply_adjoint_pde_hessian(std::vector<Real> &ahwv,
494  const std::vector<Real> &w,
495  const std::vector<Real> &v,
496  const std::vector<Real> &u,
497  const std::vector<Real> &z) const {
498  for (int i=0; i<nx_; i++) {
499  // Contribution from nonlinear term
500  ahwv[i] = 0.0;
501  if (i<nx_-1){
502  ahwv[i] += (w[i]*v[i+1] - w[i+1]*(2.0*v[i]+v[i+1]))/6.0;
503  }
504  if (i>0) {
505  ahwv[i] += (w[i-1]*(v[i-1]+2.0*v[i]) - w[i]*v[i-1])/6.0;
506  }
507  }
508  //ahwv.assign(u.size(),0.0);
509  }
510  void apply_adjoint_pde_control_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(u.size(),0.0);
516  }
517  void apply_adjoint_control_pde_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  void apply_adjoint_control_hessian(std::vector<Real> &ahwv,
525  const std::vector<Real> &w,
526  const std::vector<Real> &v,
527  const std::vector<Real> &u,
528  const std::vector<Real> &z) {
529  ahwv.assign(z.size(),0.0);
530  }
531 };
532 
533 template<class Real>
534 class L2VectorPrimal : public ROL::Vector<Real> {
535 private:
536  ROL::Ptr<std::vector<Real> > vec_;
537  ROL::Ptr<BurgersFEM<Real> > fem_;
538 
539  mutable ROL::Ptr<L2VectorDual<Real> > dual_vec_;
540 
541 public:
542  L2VectorPrimal(const ROL::Ptr<std::vector<Real> > & vec,
543  const ROL::Ptr<BurgersFEM<Real> > &fem)
544  : vec_(vec), fem_(fem), dual_vec_(ROL::nullPtr) {}
545 
546  void set( const ROL::Vector<Real> &x ) {
547  const L2VectorPrimal &ex = dynamic_cast<const L2VectorPrimal&>(x);
548  const std::vector<Real>& xval = *ex.getVector();
549  std::copy(xval.begin(),xval.end(),vec_->begin());
550  }
551 
552  void plus( const ROL::Vector<Real> &x ) {
553  const L2VectorPrimal &ex = dynamic_cast<const L2VectorPrimal&>(x);
554  const std::vector<Real>& xval = *ex.getVector();
555  unsigned dimension = vec_->size();
556  for (unsigned i=0; i<dimension; i++) {
557  (*vec_)[i] += xval[i];
558  }
559  }
560 
561  void scale( const Real alpha ) {
562  unsigned dimension = vec_->size();
563  for (unsigned i=0; i<dimension; i++) {
564  (*vec_)[i] *= alpha;
565  }
566  }
567 
568  Real dot( const ROL::Vector<Real> &x ) const {
569  const L2VectorPrimal & ex = dynamic_cast<const L2VectorPrimal&>(x);
570  const std::vector<Real>& xval = *ex.getVector();
571  return fem_->compute_L2_dot(xval,*vec_);
572  }
573 
574  Real norm() const {
575  Real val = 0;
576  val = std::sqrt( dot(*this) );
577  return val;
578  }
579 
580  ROL::Ptr<ROL::Vector<Real> > clone() const {
581  return ROL::makePtr<L2VectorPrimal>( ROL::makePtr<std::vector<Real>>(vec_->size(),0.0),fem_);
582  }
583 
584  ROL::Ptr<const std::vector<Real> > getVector() const {
585  return vec_;
586  }
587 
588  ROL::Ptr<std::vector<Real> > getVector() {
589  return vec_;
590  }
591 
592  ROL::Ptr<ROL::Vector<Real> > basis( const int i ) const {
593  ROL::Ptr<L2VectorPrimal> e
594  = ROL::makePtr<L2VectorPrimal>( ROL::makePtr<std::vector<Real>>(vec_->size(),0.0),fem_);
595  (*e->getVector())[i] = 1.0;
596  return e;
597  }
598 
599  int dimension() const {
600  return vec_->size();
601  }
602 
603  const ROL::Vector<Real>& dual() const {
604  dual_vec_ = ROL::makePtr<L2VectorDual<Real>>(
605  ROL::makePtr<std::vector<Real>>(*vec_),fem_);
606 
607  fem_->apply_mass(*(ROL::constPtrCast<std::vector<Real> >(dual_vec_->getVector())),*vec_);
608  return *dual_vec_;
609  }
610 
611 };
612 
613 template<class Real>
614 class L2VectorDual : public ROL::Vector<Real> {
615 private:
616  ROL::Ptr<std::vector<Real> > vec_;
617  ROL::Ptr<BurgersFEM<Real> > fem_;
618 
619  mutable ROL::Ptr<L2VectorPrimal<Real> > dual_vec_;
620 
621 public:
622  L2VectorDual(const ROL::Ptr<std::vector<Real> > & vec,
623  const ROL::Ptr<BurgersFEM<Real> > &fem)
624  : vec_(vec), fem_(fem), dual_vec_(ROL::nullPtr) {}
625 
626  void set( const ROL::Vector<Real> &x ) {
627  const L2VectorDual &ex = dynamic_cast<const L2VectorDual&>(x);
628  const std::vector<Real>& xval = *ex.getVector();
629  std::copy(xval.begin(),xval.end(),vec_->begin());
630  }
631 
632  void plus( const ROL::Vector<Real> &x ) {
633  const L2VectorDual &ex = dynamic_cast<const L2VectorDual&>(x);
634  const std::vector<Real>& xval = *ex.getVector();
635  unsigned dimension = vec_->size();
636  for (unsigned i=0; i<dimension; i++) {
637  (*vec_)[i] += xval[i];
638  }
639  }
640 
641  void scale( const Real alpha ) {
642  unsigned dimension = vec_->size();
643  for (unsigned i=0; i<dimension; i++) {
644  (*vec_)[i] *= alpha;
645  }
646  }
647 
648  Real dot( const ROL::Vector<Real> &x ) const {
649  const L2VectorDual & ex = dynamic_cast<const L2VectorDual&>(x);
650  const std::vector<Real>& xval = *ex.getVector();
651  unsigned dimension = vec_->size();
652  std::vector<Real> Mx(dimension,0.0);
653  fem_->apply_inverse_mass(Mx,xval);
654  Real val = 0.0;
655  for (unsigned i = 0; i < dimension; i++) {
656  val += Mx[i]*(*vec_)[i];
657  }
658  return val;
659  }
660 
661  Real norm() const {
662  Real val = 0;
663  val = std::sqrt( dot(*this) );
664  return val;
665  }
666 
667  ROL::Ptr<ROL::Vector<Real> > clone() const {
668  return ROL::makePtr<L2VectorDual>( ROL::makePtr<std::vector<Real>>(vec_->size(),0.0),fem_);
669  }
670 
671  ROL::Ptr<const std::vector<Real> > getVector() const {
672  return vec_;
673  }
674 
675  ROL::Ptr<std::vector<Real> > getVector() {
676  return vec_;
677  }
678 
679  ROL::Ptr<ROL::Vector<Real> > basis( const int i ) const {
680  ROL::Ptr<L2VectorDual> e
681  = ROL::makePtr<L2VectorDual>( ROL::makePtr<std::vector<Real>>(vec_->size(),0.0),fem_);
682  (*e->getVector())[i] = 1.0;
683  return e;
684  }
685 
686  int dimension() const {
687  return vec_->size();
688  }
689 
690  const ROL::Vector<Real>& dual() const {
691  dual_vec_ = ROL::makePtr<L2VectorPrimal<Real>>(
692  ROL::makePtr<std::vector<Real>>(*vec_),fem_);
693 
694  fem_->apply_inverse_mass(*(ROL::constPtrCast<std::vector<Real> >(dual_vec_->getVector())),*vec_);
695  return *dual_vec_;
696  }
697 
698 };
699 
700 template<class Real>
701 class H1VectorPrimal : public ROL::Vector<Real> {
702 private:
703  ROL::Ptr<std::vector<Real> > vec_;
704  ROL::Ptr<BurgersFEM<Real> > fem_;
705 
706  mutable ROL::Ptr<H1VectorDual<Real> > dual_vec_;
707 
708 public:
709  H1VectorPrimal(const ROL::Ptr<std::vector<Real> > & vec,
710  const ROL::Ptr<BurgersFEM<Real> > &fem)
711  : vec_(vec), fem_(fem), dual_vec_(ROL::nullPtr) {}
712 
713  void set( const ROL::Vector<Real> &x ) {
714  const H1VectorPrimal &ex = dynamic_cast<const H1VectorPrimal&>(x);
715  const std::vector<Real>& xval = *ex.getVector();
716  std::copy(xval.begin(),xval.end(),vec_->begin());
717  }
718 
719  void plus( const ROL::Vector<Real> &x ) {
720  const H1VectorPrimal &ex = dynamic_cast<const H1VectorPrimal&>(x);
721  const std::vector<Real>& xval = *ex.getVector();
722  unsigned dimension = vec_->size();
723  for (unsigned i=0; i<dimension; i++) {
724  (*vec_)[i] += xval[i];
725  }
726  }
727 
728  void scale( const Real alpha ) {
729  unsigned dimension = vec_->size();
730  for (unsigned i=0; i<dimension; i++) {
731  (*vec_)[i] *= alpha;
732  }
733  }
734 
735  Real dot( const ROL::Vector<Real> &x ) const {
736  const H1VectorPrimal & ex = dynamic_cast<const H1VectorPrimal&>(x);
737  const std::vector<Real>& xval = *ex.getVector();
738  return fem_->compute_H1_dot(xval,*vec_);
739  }
740 
741  Real norm() const {
742  Real val = 0;
743  val = std::sqrt( dot(*this) );
744  return val;
745  }
746 
747  ROL::Ptr<ROL::Vector<Real> > clone() const {
748  return ROL::makePtr<H1VectorPrimal>( ROL::makePtr<std::vector<Real>>(vec_->size(),0.0),fem_);
749  }
750 
751  ROL::Ptr<const std::vector<Real> > getVector() const {
752  return vec_;
753  }
754 
755  ROL::Ptr<std::vector<Real> > getVector() {
756  return vec_;
757  }
758 
759  ROL::Ptr<ROL::Vector<Real> > basis( const int i ) const {
760  ROL::Ptr<H1VectorPrimal> e
761  = ROL::makePtr<H1VectorPrimal>( ROL::makePtr<std::vector<Real>>(vec_->size(),0.0),fem_);
762  (*e->getVector())[i] = 1.0;
763  return e;
764  }
765 
766  int dimension() const {
767  return vec_->size();
768  }
769 
770  const ROL::Vector<Real>& dual() const {
771  dual_vec_ = ROL::makePtr<H1VectorDual<Real>>(
772  ROL::makePtr<std::vector<Real>>(*vec_),fem_);
773 
774  fem_->apply_H1(*(ROL::constPtrCast<std::vector<Real> >(dual_vec_->getVector())),*vec_);
775  return *dual_vec_;
776  }
777 
778 };
779 
780 template<class Real>
781 class H1VectorDual : public ROL::Vector<Real> {
782 private:
783  ROL::Ptr<std::vector<Real> > vec_;
784  ROL::Ptr<BurgersFEM<Real> > fem_;
785 
786  mutable ROL::Ptr<H1VectorPrimal<Real> > dual_vec_;
787 
788 public:
789  H1VectorDual(const ROL::Ptr<std::vector<Real> > & vec,
790  const ROL::Ptr<BurgersFEM<Real> > &fem)
791  : vec_(vec), fem_(fem), dual_vec_(ROL::nullPtr) {}
792 
793  void set( const ROL::Vector<Real> &x ) {
794  const H1VectorDual &ex = dynamic_cast<const H1VectorDual&>(x);
795  const std::vector<Real>& xval = *ex.getVector();
796  std::copy(xval.begin(),xval.end(),vec_->begin());
797  }
798 
799  void plus( const ROL::Vector<Real> &x ) {
800  const H1VectorDual &ex = dynamic_cast<const H1VectorDual&>(x);
801  const std::vector<Real>& xval = *ex.getVector();
802  unsigned dimension = vec_->size();
803  for (unsigned i=0; i<dimension; i++) {
804  (*vec_)[i] += xval[i];
805  }
806  }
807 
808  void scale( const Real alpha ) {
809  unsigned dimension = vec_->size();
810  for (unsigned i=0; i<dimension; i++) {
811  (*vec_)[i] *= alpha;
812  }
813  }
814 
815  Real dot( const ROL::Vector<Real> &x ) const {
816  const H1VectorDual & ex = dynamic_cast<const H1VectorDual&>(x);
817  const std::vector<Real>& xval = *ex.getVector();
818  unsigned dimension = vec_->size();
819  std::vector<Real> Mx(dimension,0.0);
820  fem_->apply_inverse_H1(Mx,xval);
821  Real val = 0.0;
822  for (unsigned i = 0; i < dimension; i++) {
823  val += Mx[i]*(*vec_)[i];
824  }
825  return val;
826  }
827 
828  Real norm() const {
829  Real val = 0;
830  val = std::sqrt( dot(*this) );
831  return val;
832  }
833 
834  ROL::Ptr<ROL::Vector<Real>> clone() const {
835  return ROL::makePtr<H1VectorDual>( ROL::makePtr<std::vector<Real>>(vec_->size(),0.0),fem_);
836  }
837 
838  ROL::Ptr<const std::vector<Real> > getVector() const {
839  return vec_;
840  }
841 
842  ROL::Ptr<std::vector<Real> > getVector() {
843  return vec_;
844  }
845 
846  ROL::Ptr<ROL::Vector<Real> > basis( const int i ) const {
847  ROL::Ptr<H1VectorDual> e
848  = ROL::makePtr<H1VectorDual>( ROL::makePtr<std::vector<Real>>(vec_->size(),0.0),fem_);
849  (*e->getVector())[i] = 1.0;
850  return e;
851  }
852 
853  int dimension() const {
854  return vec_->size();
855  }
856 
857  const ROL::Vector<Real>& dual() const {
858  dual_vec_ = ROL::makePtr<H1VectorPrimal<Real>>(
859  ROL::makePtr<std::vector<Real>>(*vec_),fem_);
860 
861  fem_->apply_inverse_H1(*(ROL::constPtrCast<std::vector<Real> >(dual_vec_->getVector())),*vec_);
862  return *dual_vec_;
863  }
864 
865 };
866 
867 template<class Real>
868 class Constraint_BurgersControl : public ROL::Constraint_SimOpt<Real> {
869 private:
870 
871  typedef H1VectorPrimal<Real> PrimalStateVector;
872  typedef H1VectorDual<Real> DualStateVector;
873 
874  typedef L2VectorPrimal<Real> PrimalControlVector;
875  typedef L2VectorDual<Real> DualControlVector;
876 
877  typedef H1VectorDual<Real> PrimalConstraintVector;
878  typedef H1VectorPrimal<Real> DualConstraintVector;
879 
880  ROL::Ptr<BurgersFEM<Real> > fem_;
881  bool useHessian_;
882 
883 public:
884  Constraint_BurgersControl(const ROL::Ptr<BurgersFEM<Real> > &fem,
885  const bool useHessian = true)
886  : fem_(fem), useHessian_(useHessian) {}
887 
888  void value(ROL::Vector<Real> &c, const ROL::Vector<Real> &u,
889  const ROL::Vector<Real> &z, Real &tol) {
890  ROL::Ptr<std::vector<Real> > cp =
891  dynamic_cast<PrimalConstraintVector&>(c).getVector();
892  ROL::Ptr<const std::vector<Real> > up =
893  dynamic_cast<const PrimalStateVector&>(u).getVector();
894  ROL::Ptr<const std::vector<Real> > zp =
895  dynamic_cast<const PrimalControlVector&>(z).getVector();
896 
897  fem_->compute_residual(*cp,*up,*zp);
898  }
899 
900  void applyJacobian_1(ROL::Vector<Real> &jv, const ROL::Vector<Real> &v, const ROL::Vector<Real> &u,
901  const ROL::Vector<Real> &z, Real &tol) {
902  ROL::Ptr<std::vector<Real> > jvp =
903  dynamic_cast<PrimalConstraintVector&>(jv).getVector();
904  ROL::Ptr<const std::vector<Real> > vp =
905  dynamic_cast<const PrimalStateVector&>(v).getVector();
906  ROL::Ptr<const std::vector<Real> > up =
907  dynamic_cast<const PrimalStateVector&>(u).getVector();
908  ROL::Ptr<const std::vector<Real> > zp =
909  dynamic_cast<const PrimalControlVector&>(z).getVector();
910 
911  fem_->apply_pde_jacobian(*jvp,*vp,*up,*zp);
912  }
913 
914  void applyJacobian_2(ROL::Vector<Real> &jv, const ROL::Vector<Real> &v, const ROL::Vector<Real> &u,
915  const ROL::Vector<Real> &z, Real &tol) {
916  ROL::Ptr<std::vector<Real> > jvp =
917  dynamic_cast<PrimalConstraintVector&>(jv).getVector();
918  ROL::Ptr<const std::vector<Real> > vp =
919  dynamic_cast<const PrimalControlVector&>(v).getVector();
920  ROL::Ptr<const std::vector<Real> > up =
921  dynamic_cast<const PrimalStateVector&>(u).getVector();
922  ROL::Ptr<const std::vector<Real> > zp =
923  dynamic_cast<const PrimalControlVector&>(z).getVector();
924 
925  fem_->apply_control_jacobian(*jvp,*vp,*up,*zp);
926  }
927 
928  void applyInverseJacobian_1(ROL::Vector<Real> &ijv, const ROL::Vector<Real> &v, const ROL::Vector<Real> &u,
929  const ROL::Vector<Real> &z, Real &tol) {
930  ROL::Ptr<std::vector<Real> > ijvp =
931  dynamic_cast<PrimalStateVector&>(ijv).getVector();
932  ROL::Ptr<const std::vector<Real> > vp =
933  dynamic_cast<const PrimalConstraintVector&>(v).getVector();
934  ROL::Ptr<const std::vector<Real> > up =
935  dynamic_cast<const PrimalStateVector&>(u).getVector();
936  ROL::Ptr<const std::vector<Real> > zp =
937  dynamic_cast<const PrimalControlVector&>(z).getVector();
938 
939  fem_->apply_inverse_pde_jacobian(*ijvp,*vp,*up,*zp);
940  }
941 
942  void applyAdjointJacobian_1(ROL::Vector<Real> &ajv, const ROL::Vector<Real> &v, const ROL::Vector<Real> &u,
943  const ROL::Vector<Real> &z, Real &tol) {
944  ROL::Ptr<std::vector<Real> > jvp =
945  dynamic_cast<DualStateVector&>(ajv).getVector();
946  ROL::Ptr<const std::vector<Real> > vp =
947  dynamic_cast<const DualConstraintVector&>(v).getVector();
948  ROL::Ptr<const std::vector<Real> > up =
949  dynamic_cast<const PrimalStateVector&>(u).getVector();
950  ROL::Ptr<const std::vector<Real> > zp =
951  dynamic_cast<const PrimalControlVector&>(z).getVector();
952 
953  fem_->apply_adjoint_pde_jacobian(*jvp,*vp,*up,*zp);
954  }
955 
956  void applyAdjointJacobian_2(ROL::Vector<Real> &jv, const ROL::Vector<Real> &v, const ROL::Vector<Real> &u,
957  const ROL::Vector<Real> &z, Real &tol) {
958  ROL::Ptr<std::vector<Real> > jvp =
959  dynamic_cast<DualControlVector&>(jv).getVector();
960  ROL::Ptr<const std::vector<Real> > vp =
961  dynamic_cast<const DualConstraintVector&>(v).getVector();
962  ROL::Ptr<const std::vector<Real> > up =
963  dynamic_cast<const PrimalStateVector&>(u).getVector();
964  ROL::Ptr<const std::vector<Real> > zp =
965  dynamic_cast<const PrimalControlVector&>(z).getVector();
966 
967  fem_->apply_adjoint_control_jacobian(*jvp,*vp,*up,*zp);
968  }
969 
970  void applyInverseAdjointJacobian_1(ROL::Vector<Real> &iajv, const ROL::Vector<Real> &v,
971  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
972  ROL::Ptr<std::vector<Real> > iajvp =
973  dynamic_cast<DualConstraintVector&>(iajv).getVector();
974  ROL::Ptr<const std::vector<Real> > vp =
975  dynamic_cast<const DualStateVector&>(v).getVector();
976  ROL::Ptr<const std::vector<Real> > up =
977  dynamic_cast<const PrimalStateVector&>(u).getVector();
978  ROL::Ptr<const std::vector<Real> > zp =
979  dynamic_cast<const PrimalControlVector&>(z).getVector();
980 
981  fem_->apply_inverse_adjoint_pde_jacobian(*iajvp,*vp,*up,*zp);
982  }
983 
984  void applyAdjointHessian_11(ROL::Vector<Real> &ahwv, const ROL::Vector<Real> &w, const ROL::Vector<Real> &v,
985  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
986  if ( useHessian_ ) {
987  ROL::Ptr<std::vector<Real> > ahwvp =
988  dynamic_cast<DualStateVector&>(ahwv).getVector();
989  ROL::Ptr<const std::vector<Real> > wp =
990  dynamic_cast<const DualConstraintVector&>(w).getVector();
991  ROL::Ptr<const std::vector<Real> > vp =
992  dynamic_cast<const PrimalStateVector&>(v).getVector();
993  ROL::Ptr<const std::vector<Real> > up =
994  dynamic_cast<const PrimalStateVector&>(u).getVector();
995  ROL::Ptr<const std::vector<Real> > zp =
996  dynamic_cast<const PrimalControlVector&>(z).getVector();
997 
998  fem_->apply_adjoint_pde_hessian(*ahwvp,*wp,*vp,*up,*zp);
999  }
1000  else {
1001  ahwv.zero();
1002  }
1003  }
1004 
1005  void applyAdjointHessian_12(ROL::Vector<Real> &ahwv, const ROL::Vector<Real> &w, const ROL::Vector<Real> &v,
1006  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
1007  if ( useHessian_ ) {
1008  ROL::Ptr<std::vector<Real> > ahwvp =
1009  dynamic_cast<DualControlVector&>(ahwv).getVector();
1010  ROL::Ptr<const std::vector<Real> > wp =
1011  dynamic_cast<const DualConstraintVector&>(w).getVector();
1012  ROL::Ptr<const std::vector<Real> > vp =
1013  dynamic_cast<const PrimalStateVector&>(v).getVector();
1014  ROL::Ptr<const std::vector<Real> > up =
1015  dynamic_cast<const PrimalStateVector&>(u).getVector();
1016  ROL::Ptr<const std::vector<Real> > zp =
1017  dynamic_cast<const PrimalControlVector&>(z).getVector();
1018 
1019  fem_->apply_adjoint_control_pde_hessian(*ahwvp,*wp,*vp,*up,*zp);
1020  }
1021  else {
1022  ahwv.zero();
1023  }
1024  }
1025  void applyAdjointHessian_21(ROL::Vector<Real> &ahwv, const ROL::Vector<Real> &w, const ROL::Vector<Real> &v,
1026  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
1027  if ( useHessian_ ) {
1028  ROL::Ptr<std::vector<Real> > ahwvp =
1029  dynamic_cast<DualStateVector&>(ahwv).getVector();
1030  ROL::Ptr<const std::vector<Real> > wp =
1031  dynamic_cast<const DualConstraintVector&>(w).getVector();
1032  ROL::Ptr<const std::vector<Real> > vp =
1033  dynamic_cast<const PrimalControlVector&>(v).getVector();
1034  ROL::Ptr<const std::vector<Real> > up =
1035  dynamic_cast<const PrimalStateVector&>(u).getVector();
1036  ROL::Ptr<const std::vector<Real> > zp =
1037  dynamic_cast<const PrimalControlVector&>(z).getVector();
1038 
1039  fem_->apply_adjoint_pde_control_hessian(*ahwvp,*wp,*vp,*up,*zp);
1040  }
1041  else {
1042  ahwv.zero();
1043  }
1044  }
1045  void applyAdjointHessian_22(ROL::Vector<Real> &ahwv, const ROL::Vector<Real> &w, const ROL::Vector<Real> &v,
1046  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
1047  if ( useHessian_ ) {
1048  ROL::Ptr<std::vector<Real> > ahwvp =
1049  dynamic_cast<DualControlVector&>(ahwv).getVector();
1050  ROL::Ptr<const std::vector<Real> > wp =
1051  dynamic_cast<const DualConstraintVector&>(w).getVector();
1052  ROL::Ptr<const std::vector<Real> > vp =
1053  dynamic_cast<const PrimalControlVector&>(v).getVector();
1054  ROL::Ptr<const std::vector<Real> > up =
1055  dynamic_cast<const PrimalStateVector&>(u).getVector();
1056  ROL::Ptr<const std::vector<Real> > zp =
1057  dynamic_cast<const PrimalControlVector&>(z).getVector();
1058 
1059  fem_->apply_adjoint_control_hessian(*ahwvp,*wp,*vp,*up,*zp);
1060  }
1061  else {
1062  ahwv.zero();
1063  }
1064  }
1065 };
BurgersFEM::compute_H1_norm
Real compute_H1_norm(const std::vector< Real > &r) const
Definition: test_04.hpp:268
L2VectorPrimal
Definition: test_04.hpp:57
L2VectorPrimal::scale
void scale(const Real alpha)
Compute where .
Definition: test_04.hpp:561
H1VectorDual::basis
ROL::Ptr< ROL::Vector< Real > > basis(const int i) const
Return i-th basis vector.
Definition: test_04.hpp:846
H1VectorDual::clone
ROL::Ptr< ROL::Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
Definition: test_04.hpp:834
L2VectorDual::scale
void scale(const Real alpha)
Compute where .
Definition: test_04.hpp:641
BurgersFEM::mesh_spacing
Real mesh_spacing(void) const
Definition: test_04.hpp:150
L2VectorDual::L2VectorDual
L2VectorDual(const ROL::Ptr< std::vector< Real > > &vec, const ROL::Ptr< BurgersFEM< Real > > &fem)
Definition: test_04.hpp:622
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: test_04.hpp:418
H1VectorDual::plus
void plus(const ROL::Vector< Real > &x)
Compute , where .
Definition: test_04.hpp:799
H1VectorDual
Definition: test_04.hpp:66
ROL_Objective_SimOpt.hpp
H1VectorPrimal::dual_vec_
ROL::Ptr< H1VectorDual< Real > > dual_vec_
Definition: test_04.hpp:706
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: test_04.hpp:524
BurgersFEM
Definition: test_04.hpp:69
L2VectorDual::getVector
ROL::Ptr< std::vector< Real > > getVector()
Definition: test_04.hpp:675
L2VectorDual::vec_
ROL::Ptr< std::vector< Real > > vec_
Definition: test_04.hpp:616
L2VectorDual::plus
void plus(const ROL::Vector< Real > &x)
Compute , where .
Definition: test_04.hpp:632
L2VectorPrimal::clone
ROL::Ptr< ROL::Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
Definition: test_04.hpp:580
L2VectorPrimal::L2VectorPrimal
L2VectorPrimal(const ROL::Ptr< std::vector< Real > > &vec, const ROL::Ptr< BurgersFEM< Real > > &fem)
Definition: test_04.hpp:542
BurgersFEM::set_problem_data
void set_problem_data(const Real nu, const Real f, const Real u0, const Real u1)
Definition: test_04.hpp:139
L2VectorPrimal::basis
ROL::Ptr< ROL::Vector< Real > > basis(const int i) const
Return i-th basis vector.
Definition: test_04.hpp:592
H1VectorPrimal::scale
void scale(const Real alpha)
Compute where .
Definition: test_04.hpp:728
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: test_04.hpp:100
L2VectorPrimal::vec_
ROL::Ptr< std::vector< Real > > vec_
Definition: test_04.hpp:536
L2VectorPrimal::getVector
ROL::Ptr< const std::vector< Real > > getVector() const
Definition: test_04.hpp:584
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: test_04.hpp:928
L2VectorDual::fem_
ROL::Ptr< BurgersFEM< Real > > fem_
Definition: test_04.hpp:617
H1VectorDual::dot
Real dot(const ROL::Vector< Real > &x) const
Compute where .
Definition: test_04.hpp:815
BurgersFEM::nu_
Real nu_
Definition: test_04.hpp:73
ROL::Vector::zero
virtual void zero()
Set to zero vector.
Definition: ROL_Vector.hpp:167
Constraint_BurgersControl::PrimalStateVector
H1VectorPrimal< Real > PrimalStateVector
Definition: test_04.hpp:871
Constraint_BurgersControl
Definition: test_04.hpp:868
BurgersFEM::apply_inverse_mass
void apply_inverse_mass(std::vector< Real > &Mu, const std::vector< Real > &u) const
Definition: test_04.hpp:198
BurgersFEM::test_inverse_H1
void test_inverse_H1(std::ostream &outStream=std::cout)
Definition: test_04.hpp:300
Constraint_BurgersControl::PrimalConstraintVector
H1VectorDual< Real > PrimalConstraintVector
Definition: test_04.hpp:877
H1VectorDual::getVector
ROL::Ptr< std::vector< Real > > getVector()
Definition: test_04.hpp:842
L2VectorPrimal::fem_
ROL::Ptr< BurgersFEM< Real > > fem_
Definition: test_04.hpp:537
BurgersFEM::update
void update(std::vector< Real > &u, const std::vector< Real > &s, const Real alpha=1.0) const
Definition: test_04.hpp:82
H1VectorPrimal::fem_
ROL::Ptr< BurgersFEM< Real > > fem_
Definition: test_04.hpp:704
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: test_04.hpp:439
H1VectorPrimal::H1VectorPrimal
H1VectorPrimal(const ROL::Ptr< std::vector< Real > > &vec, const ROL::Ptr< BurgersFEM< Real > > &fem)
Definition: test_04.hpp:709
H1VectorPrimal::dot
Real dot(const ROL::Vector< Real > &x) const
Compute where .
Definition: test_04.hpp:735
L2VectorPrimal::dot
Real dot(const ROL::Vector< Real > &x) const
Compute where .
Definition: test_04.hpp:568
BurgersFEM::f_
Real f_
Definition: test_04.hpp:77
L2VectorDual::set
void set(const ROL::Vector< Real > &x)
Set where .
Definition: test_04.hpp:626
BurgersFEM::compute_L2_norm
Real compute_L2_norm(const std::vector< Real > &r) const
Definition: test_04.hpp:176
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: test_04.hpp:857
H1VectorDual::dimension
int dimension() const
Return dimension of the vector space.
Definition: test_04.hpp:853
H1VectorPrimal::dimension
int dimension() const
Return dimension of the vector space.
Definition: test_04.hpp:766
H1VectorPrimal::norm
Real norm() const
Returns where .
Definition: test_04.hpp:741
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: test_04.hpp:517
H1VectorPrimal::getVector
ROL::Ptr< const std::vector< Real > > getVector() const
Definition: test_04.hpp:751
BurgersFEM::num_dof
int num_dof(void) const
Definition: test_04.hpp:146
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: test_04.hpp:942
H1VectorPrimal
Definition: test_04.hpp:63
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: test_04.hpp:603
L2VectorPrimal::plus
void plus(const ROL::Vector< Real > &x)
Compute , where .
Definition: test_04.hpp:552
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: test_04.hpp:770
BurgersFEM::BurgersFEM
BurgersFEM(int nx=128, Real nl=1.0, Real cH1=1.0, Real cL2=1.0)
Definition: test_04.hpp:131
L2VectorDual
Definition: test_04.hpp:60
ROL_Constraint_SimOpt.hpp
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: test_04.hpp:385
BurgersFEM::cL2_
Real cL2_
Definition: test_04.hpp:79
H1VectorPrimal::getVector
ROL::Ptr< std::vector< Real > > getVector()
Definition: test_04.hpp:755
H1VectorDual::set
void set(const ROL::Vector< Real > &x)
Set where .
Definition: test_04.hpp:793
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: test_04.hpp:984
ROL::Vector
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:80
H1VectorDual::dual_vec_
ROL::Ptr< H1VectorPrimal< Real > > dual_vec_
Definition: test_04.hpp:786
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: test_04.hpp:510
Constraint_BurgersControl::useHessian_
bool useHessian_
Definition: test_04.hpp:881
BurgersFEM::dx_
Real dx_
Definition: test_04.hpp:72
H1VectorDual::norm
Real norm() const
Returns where .
Definition: test_04.hpp:828
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: test_04.hpp:1005
L2VectorPrimal::norm
Real norm() const
Returns where .
Definition: test_04.hpp:574
L2VectorDual::dual_vec_
ROL::Ptr< L2VectorPrimal< Real > > dual_vec_
Definition: test_04.hpp:619
H1VectorDual::H1VectorDual
H1VectorDual(const ROL::Ptr< std::vector< Real > > &vec, const ROL::Ptr< BurgersFEM< Real > > &fem)
Definition: test_04.hpp:789
H1VectorPrimal::plus
void plus(const ROL::Vector< Real > &x)
Compute , where .
Definition: test_04.hpp:719
ROL_Types.hpp
Contains definitions of custom data types in ROL.
BurgersFEM::compute_residual
void compute_residual(std::vector< Real > &r, const std::vector< Real > &u, const std::vector< Real > &z) const
Definition: test_04.hpp:321
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: test_04.hpp:900
BurgersFEM::u0_
Real u0_
Definition: test_04.hpp:75
L2VectorDual::norm
Real norm() const
Returns where .
Definition: test_04.hpp:661
H1VectorDual::getVector
ROL::Ptr< const std::vector< Real > > getVector() const
Definition: test_04.hpp:838
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: test_04.hpp:1025
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: test_04.hpp:1045
BurgersFEM::apply_inverse_H1
void apply_inverse_H1(std::vector< Real > &Mu, const std::vector< Real > &u) const
Definition: test_04.hpp:292
BurgersFEM::axpy
void axpy(std::vector< Real > &out, const Real a, const std::vector< Real > &x, const std::vector< Real > &y) const
Definition: test_04.hpp:88
Constraint_BurgersControl::Constraint_BurgersControl
Constraint_BurgersControl(const ROL::Ptr< BurgersFEM< Real > > &fem, const bool useHessian=true)
Definition: test_04.hpp:884
ROL
Definition: ROL_ElementwiseVector.hpp:61
H1VectorDual::fem_
ROL::Ptr< BurgersFEM< Real > > fem_
Definition: test_04.hpp:784
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: test_04.hpp:357
ROL::Constraint_SimOpt
Defines the constraint operator interface for simulation-based optimization.
Definition: ROL_Constraint_SimOpt.hpp:100
Constraint_BurgersControl::DualStateVector
H1VectorDual< Real > DualStateVector
Definition: test_04.hpp:872
H1VectorPrimal::vec_
ROL::Ptr< std::vector< Real > > vec_
Definition: test_04.hpp:703
Constraint_BurgersControl::DualConstraintVector
H1VectorPrimal< Real > DualConstraintVector
Definition: test_04.hpp:878
BurgersFEM::nx_
int nx_
Definition: test_04.hpp:71
L2VectorDual::dot
Real dot(const ROL::Vector< Real > &x) const
Compute where .
Definition: test_04.hpp:648
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: test_04.hpp:970
L2VectorPrimal::getVector
ROL::Ptr< std::vector< Real > > getVector()
Definition: test_04.hpp:588
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: test_04.hpp:914
L2VectorDual::basis
ROL::Ptr< ROL::Vector< Real > > basis(const int i) const
Return i-th basis vector.
Definition: test_04.hpp:679
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: test_04.hpp:404
L2VectorPrimal::set
void set(const ROL::Vector< Real > &x)
Set where .
Definition: test_04.hpp:546
H1VectorDual::scale
void scale(const Real alpha)
Compute where .
Definition: test_04.hpp:808
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: test_04.hpp:456
BurgersFEM::scale
void scale(std::vector< Real > &u, const Real alpha=0.0) const
Definition: test_04.hpp:94
BurgersFEM::compute_H1_dot
Real compute_H1_dot(const std::vector< Real > &x, const std::vector< Real > &y) const
Definition: test_04.hpp:248
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: test_04.hpp:956
BurgersFEM::apply_H1
void apply_H1(std::vector< Real > &Mu, const std::vector< Real > &u) const
Definition: test_04.hpp:273
BurgersFEM::nl_
Real nl_
Definition: test_04.hpp:74
Constraint_BurgersControl::DualControlVector
L2VectorDual< Real > DualControlVector
Definition: test_04.hpp:875
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: test_04.hpp:888
BurgersFEM::u1_
Real u1_
Definition: test_04.hpp:76
BurgersFEM::compute_L2_dot
Real compute_L2_dot(const std::vector< Real > &x, const std::vector< Real > &y) const
Definition: test_04.hpp:158
H1VectorPrimal::set
void set(const ROL::Vector< Real > &x)
Set where .
Definition: test_04.hpp:713
ROL_Vector.hpp
H1VectorPrimal::basis
ROL::Ptr< ROL::Vector< Real > > basis(const int i) const
Return i-th basis vector.
Definition: test_04.hpp:759
BurgersFEM::test_inverse_mass
void test_inverse_mass(std::ostream &outStream=std::cout)
Definition: test_04.hpp:211
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: test_04.hpp:493
L2VectorDual::dimension
int dimension() const
Return dimension of the vector space.
Definition: test_04.hpp:686
H1VectorPrimal::clone
ROL::Ptr< ROL::Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
Definition: test_04.hpp:747
H1VectorDual::vec_
ROL::Ptr< std::vector< Real > > vec_
Definition: test_04.hpp:783
BurgersFEM::apply_mass
void apply_mass(std::vector< Real > &Mu, const std::vector< Real > &u) const
Definition: test_04.hpp:181
L2VectorPrimal::dual_vec_
ROL::Ptr< L2VectorDual< Real > > dual_vec_
Definition: test_04.hpp:539
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: test_04.hpp:467
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: test_04.hpp:690
L2VectorDual::clone
ROL::Ptr< ROL::Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
Definition: test_04.hpp:667
Constraint_BurgersControl::PrimalControlVector
L2VectorPrimal< Real > PrimalControlVector
Definition: test_04.hpp:874
L2VectorPrimal::dimension
int dimension() const
Return dimension of the vector space.
Definition: test_04.hpp:599
ROL_BoundConstraint.hpp
BurgersFEM::cH1_
Real cH1_
Definition: test_04.hpp:78
Constraint_BurgersControl::fem_
ROL::Ptr< BurgersFEM< Real > > fem_
Definition: test_04.hpp:880
L2VectorDual::getVector
ROL::Ptr< const std::vector< Real > > getVector() const
Definition: test_04.hpp:671
Generated by   doxygen 1.8.16