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