ROL
test_04.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Rapid Optimization Library (ROL) Package
4 //
5 // Copyright 2014 NTESS and the ROL contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
14 #include "ROL_Types.hpp"
15 #include "ROL_StdVector.hpp"
16 #include "ROL_BoundConstraint.hpp"
18 #include "ROL_Objective_SimOpt.hpp"
19 
20 #include "Teuchos_LAPACK.hpp"
21 
22 template<class Real>
24 
25 template<class Real>
27 
28 template<class Real>
30 
31 template<class Real>
33 
34 template<class Real>
35 class BurgersFEM {
36 private:
37  int nx_;
38  Real dx_;
39  Real nu_;
40  Real nl_;
41  Real u0_;
42  Real u1_;
43  Real f_;
44  Real cH1_;
45  Real cL2_;
46 
47 private:
48  void update(std::vector<Real> &u, const std::vector<Real> &s, const Real alpha=1.0) const {
49  for (unsigned i=0; i<u.size(); i++) {
50  u[i] += alpha*s[i];
51  }
52  }
53 
54  void axpy(std::vector<Real> &out, const Real a, const std::vector<Real> &x, const std::vector<Real> &y) const {
55  for (unsigned i=0; i < x.size(); i++) {
56  out[i] = a*x[i] + y[i];
57  }
58  }
59 
60  void scale(std::vector<Real> &u, const Real alpha=0.0) const {
61  for (unsigned i=0; i<u.size(); i++) {
62  u[i] *= alpha;
63  }
64  }
65 
66  void linear_solve(std::vector<Real> &u, std::vector<Real> &dl, std::vector<Real> &d, std::vector<Real> &du,
67  const std::vector<Real> &r, const bool transpose = false) const {
68  if ( r.size() == 1 ) {
69  u.resize(1,r[0]/d[0]);
70  }
71  else if ( r.size() == 2 ) {
72  u.resize(2,0.0);
73  Real det = d[0]*d[1] - dl[0]*du[0];
74  u[0] = (d[1]*r[0] - du[0]*r[1])/det;
75  u[1] = (d[0]*r[1] - dl[0]*r[0])/det;
76  }
77  else {
78  u.assign(r.begin(),r.end());
79  // Perform LDL factorization
80  Teuchos::LAPACK<int,Real> lp;
81  std::vector<Real> du2(r.size()-2,0.0);
82  std::vector<int> ipiv(r.size(),0);
83  int info;
84  int dim = r.size();
85  int ldb = r.size();
86  int nhrs = 1;
87  lp.GTTRF(dim,&dl[0],&d[0],&du[0],&du2[0],&ipiv[0],&info);
88  char trans = 'N';
89  if ( transpose ) {
90  trans = 'T';
91  }
92  lp.GTTRS(trans,dim,nhrs,&dl[0],&d[0],&du[0],&du2[0],&ipiv[0],&u[0],ldb,&info);
93  }
94  }
95 
96 public:
97  BurgersFEM(int nx = 128, Real nl = 1.0, Real cH1 = 1.0, Real cL2 = 1.0)
98  : nx_(nx), dx_(1.0/((Real)nx+1.0)), nl_(nl), cH1_(cH1), cL2_(cL2) {
99  nu_ = 1.e-2;
100  f_ = 0.0;
101  u0_ = 1.0;
102  u1_ = 0.0;
103  }
104 
105  void set_problem_data(const Real nu, const Real f, const Real u0, const Real u1) {
106  nu_ = std::pow(10.0,nu-2.0);
107  f_ = 0.01*f;
108  u0_ = 1.0+0.001*u0;
109  u1_ = 0.001*u1;
110  }
111 
112  int num_dof(void) const {
113  return nx_;
114  }
115 
116  Real mesh_spacing(void) const {
117  return dx_;
118  }
119 
120  /***************************************************************************/
121  /*********************** L2 INNER PRODUCT FUNCTIONS ************************/
122  /***************************************************************************/
123  // Compute L2 inner product
124  Real compute_L2_dot(const std::vector<Real> &x, const std::vector<Real> &y) const {
125  Real ip = 0.0;
126  Real c = (((int)x.size()==nx_) ? 4.0 : 2.0);
127  for (unsigned i=0; i<x.size(); i++) {
128  if ( i == 0 ) {
129  ip += dx_/6.0*(c*x[i] + x[i+1])*y[i];
130  }
131  else if ( i == x.size()-1 ) {
132  ip += dx_/6.0*(x[i-1] + c*x[i])*y[i];
133  }
134  else {
135  ip += dx_/6.0*(x[i-1] + 4.0*x[i] + x[i+1])*y[i];
136  }
137  }
138  return ip;
139  }
140 
141  // compute L2 norm
142  Real compute_L2_norm(const std::vector<Real> &r) const {
143  return std::sqrt(compute_L2_dot(r,r));
144  }
145 
146  // Apply L2 Reisz operator
147  void apply_mass(std::vector<Real> &Mu, const std::vector<Real> &u ) const {
148  Mu.resize(u.size(),0.0);
149  Real c = (((int)u.size()==nx_) ? 4.0 : 2.0);
150  for (unsigned i=0; i<u.size(); i++) {
151  if ( i == 0 ) {
152  Mu[i] = dx_/6.0*(c*u[i] + u[i+1]);
153  }
154  else if ( i == u.size()-1 ) {
155  Mu[i] = dx_/6.0*(u[i-1] + c*u[i]);
156  }
157  else {
158  Mu[i] = dx_/6.0*(u[i-1] + 4.0*u[i] + u[i+1]);
159  }
160  }
161  }
162 
163  // Apply L2 inverse Reisz operator
164  void apply_inverse_mass(std::vector<Real> &Mu, const std::vector<Real> &u) const {
165  unsigned nx = u.size();
166  // Build mass matrix
167  std::vector<Real> dl(nx-1,dx_/6.0);
168  std::vector<Real> d(nx,2.0*dx_/3.0);
169  std::vector<Real> du(nx-1,dx_/6.0);
170  if ( (int)nx != nx_ ) {
171  d[ 0] = dx_/3.0;
172  d[nx-1] = dx_/3.0;
173  }
174  linear_solve(Mu,dl,d,du,u);
175  }
176 
177  void test_inverse_mass(std::ostream &outStream = std::cout) {
178  // State Mass Matrix
179  std::vector<Real> u(nx_,0.0), Mu(nx_,0.0), iMMu(nx_,0.0), diff(nx_,0.0);
180  for (int i = 0; i < nx_; i++) {
181  u[i] = 2.0*(Real)rand()/(Real)RAND_MAX - 1.0;
182  }
183  apply_mass(Mu,u);
184  apply_inverse_mass(iMMu,Mu);
185  axpy(diff,-1.0,iMMu,u);
186  Real error = compute_L2_norm(diff);
187  Real normu = compute_L2_norm(u);
188  outStream << "\nTest Inverse State Mass Matrix\n";
189  outStream << " ||u - inv(M)Mu|| = " << error << "\n";
190  outStream << " ||u|| = " << normu << "\n";
191  outStream << " Relative Error = " << error/normu << "\n";
192  outStream << "\n";
193  // Control Mass Matrix
194  u.resize(nx_+2,0.0); Mu.resize(nx_+2,0.0); iMMu.resize(nx_+2,0.0); diff.resize(nx_+2,0.0);
195  for (int i = 0; i < nx_+2; i++) {
196  u[i] = 2.0*(Real)rand()/(Real)RAND_MAX - 1.0;
197  }
198  apply_mass(Mu,u);
199  apply_inverse_mass(iMMu,Mu);
200  axpy(diff,-1.0,iMMu,u);
201  error = compute_L2_norm(diff);
202  normu = compute_L2_norm(u);
203  outStream << "\nTest Inverse Control Mass Matrix\n";
204  outStream << " ||z - inv(M)Mz|| = " << error << "\n";
205  outStream << " ||z|| = " << normu << "\n";
206  outStream << " Relative Error = " << error/normu << "\n";
207  outStream << "\n";
208  }
209 
210  /***************************************************************************/
211  /*********************** H1 INNER PRODUCT FUNCTIONS ************************/
212  /***************************************************************************/
213  // Compute H1 inner product
214  Real compute_H1_dot(const std::vector<Real> &x, const std::vector<Real> &y) const {
215  Real ip = 0.0;
216  for (int i=0; i<nx_; i++) {
217  if ( i == 0 ) {
218  ip += cL2_*dx_/6.0*(4.0*x[i] + x[i+1])*y[i]; // Mass term
219  ip += cH1_*(2.0*x[i] - x[i+1])/dx_*y[i]; // Stiffness term
220  }
221  else if ( i == nx_-1 ) {
222  ip += cL2_*dx_/6.0*(x[i-1] + 4.0*x[i])*y[i]; // Mass term
223  ip += cH1_*(2.0*x[i] - x[i-1])/dx_*y[i]; // Stiffness term
224  }
225  else {
226  ip += cL2_*dx_/6.0*(x[i-1] + 4.0*x[i] + x[i+1])*y[i]; // Mass term
227  ip += cH1_*(2.0*x[i] - x[i-1] - x[i+1])/dx_*y[i]; // Stiffness term
228  }
229  }
230  return ip;
231  }
232 
233  // compute H1 norm
234  Real compute_H1_norm(const std::vector<Real> &r) const {
235  return std::sqrt(compute_H1_dot(r,r));
236  }
237 
238  // Apply H2 Reisz operator
239  void apply_H1(std::vector<Real> &Mu, const std::vector<Real> &u ) const {
240  Mu.resize(nx_,0.0);
241  for (int i=0; i<nx_; i++) {
242  if ( i == 0 ) {
243  Mu[i] = cL2_*dx_/6.0*(4.0*u[i] + u[i+1])
244  + cH1_*(2.0*u[i] - u[i+1])/dx_;
245  }
246  else if ( i == nx_-1 ) {
247  Mu[i] = cL2_*dx_/6.0*(u[i-1] + 4.0*u[i])
248  + cH1_*(2.0*u[i] - u[i-1])/dx_;
249  }
250  else {
251  Mu[i] = cL2_*dx_/6.0*(u[i-1] + 4.0*u[i] + u[i+1])
252  + cH1_*(2.0*u[i] - u[i-1] - u[i+1])/dx_;
253  }
254  }
255  }
256 
257  // Apply H1 inverse Reisz operator
258  void apply_inverse_H1(std::vector<Real> &Mu, const std::vector<Real> &u) const {
259  // Build mass matrix
260  std::vector<Real> dl(nx_-1,cL2_*dx_/6.0 - cH1_/dx_);
261  std::vector<Real> d(nx_,2.0*(cL2_*dx_/3.0 + cH1_/dx_));
262  std::vector<Real> du(nx_-1,cL2_*dx_/6.0 - cH1_/dx_);
263  linear_solve(Mu,dl,d,du,u);
264  }
265 
266  void test_inverse_H1(std::ostream &outStream = std::cout) {
267  std::vector<Real> u(nx_,0.0), Mu(nx_,0.0), iMMu(nx_,0.0), diff(nx_,0.0);
268  for (int i = 0; i < nx_; i++) {
269  u[i] = 2.0*(Real)rand()/(Real)RAND_MAX - 1.0;
270  }
271  apply_H1(Mu,u);
272  apply_inverse_H1(iMMu,Mu);
273  axpy(diff,-1.0,iMMu,u);
274  Real error = compute_H1_norm(diff);
275  Real normu = compute_H1_norm(u);
276  outStream << "\nTest Inverse State H1 Matrix\n";
277  outStream << " ||u - inv(M)Mu|| = " << error << "\n";
278  outStream << " ||u|| = " << normu << "\n";
279  outStream << " Relative Error = " << error/normu << "\n";
280  outStream << "\n";
281  }
282 
283  /***************************************************************************/
284  /*********************** PDE RESIDUAL AND SOLVE ****************************/
285  /***************************************************************************/
286  // Compute current PDE residual
287  void compute_residual(std::vector<Real> &r, const std::vector<Real> &u,
288  const std::vector<Real> &z) const {
289  r.clear();
290  r.resize(nx_,0.0);
291  for (int i=0; i<nx_; i++) {
292  // Contribution from stiffness term
293  if ( i==0 ) {
294  r[i] = nu_/dx_*(2.0*u[i]-u[i+1]);
295  }
296  else if (i==nx_-1) {
297  r[i] = nu_/dx_*(2.0*u[i]-u[i-1]);
298  }
299  else {
300  r[i] = nu_/dx_*(2.0*u[i]-u[i-1]-u[i+1]);
301  }
302  // Contribution from nonlinear term
303  if (i<nx_-1){
304  r[i] += nl_*u[i+1]*(u[i]+u[i+1])/6.0;
305  }
306  if (i>0) {
307  r[i] -= nl_*u[i-1]*(u[i-1]+u[i])/6.0;
308  }
309  // Contribution from control
310  r[i] -= dx_/6.0*(z[i]+4.0*z[i+1]+z[i+2]);
311  // Contribution from right-hand side
312  r[i] -= dx_*f_;
313  }
314  // Contribution from Dirichlet boundary terms
315  r[0] -= nl_*(u0_*u[ 0]/6.0 + u0_*u0_/6.0) + nu_*u0_/dx_;
316  r[nx_-1] += nl_*(u1_*u[nx_-1]/6.0 + u1_*u1_/6.0) - nu_*u1_/dx_;
317  }
318 
319  /***************************************************************************/
320  /*********************** PDE JACOBIAN FUNCTIONS ****************************/
321  /***************************************************************************/
322  // Build PDE Jacobian trigiagonal matrix
323  void compute_pde_jacobian(std::vector<Real> &dl, // Lower diagonal
324  std::vector<Real> &d, // Diagonal
325  std::vector<Real> &du, // Upper diagonal
326  const std::vector<Real> &u) const { // State variable
327  // Get Diagonal and Off-Diagonal Entries of linear PDE Jacobian
328  d.clear();
329  d.resize(nx_,nu_*2.0/dx_);
330  dl.clear();
331  dl.resize(nx_-1,-nu_/dx_);
332  du.clear();
333  du.resize(nx_-1,-nu_/dx_);
334  // Contribution from nonlinearity
335  for (int i=0; i<nx_; i++) {
336  if (i<nx_-1) {
337  dl[i] += nl_*(-2.0*u[i]-u[i+1])/6.0;
338  d[i] += nl_*u[i+1]/6.0;
339  }
340  if (i>0) {
341  d[i] -= nl_*u[i-1]/6.0;
342  du[i-1] += nl_*(u[i-1]+2.0*u[i])/6.0;
343  }
344  }
345  // Contribution from Dirichlet boundary conditions
346  d[ 0] -= nl_*u0_/6.0;
347  d[nx_-1] += nl_*u1_/6.0;
348  }
349 
350  // Apply PDE Jacobian to a vector
351  void apply_pde_jacobian(std::vector<Real> &jv,
352  const std::vector<Real> &v,
353  const std::vector<Real> &u,
354  const std::vector<Real> &z) const {
355  // Fill jv
356  for (int i = 0; i < nx_; i++) {
357  jv[i] = nu_/dx_*2.0*v[i];
358  if ( i > 0 ) {
359  jv[i] += -nu_/dx_*v[i-1]-nl_*(u[i-1]/6.0*v[i]+(u[i]+2.0*u[i-1])/6.0*v[i-1]);
360  }
361  if ( i < nx_-1 ) {
362  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]);
363  }
364  }
365  jv[ 0] -= nl_*u0_/6.0*v[0];
366  jv[nx_-1] += nl_*u1_/6.0*v[nx_-1];
367  }
368 
369  // Apply inverse PDE Jacobian to a vector
370  void apply_inverse_pde_jacobian(std::vector<Real> &ijv,
371  const std::vector<Real> &v,
372  const std::vector<Real> &u,
373  const std::vector<Real> &z) const {
374  // Get PDE Jacobian
375  std::vector<Real> d(nx_,0.0);
376  std::vector<Real> dl(nx_-1,0.0);
377  std::vector<Real> du(nx_-1,0.0);
378  compute_pde_jacobian(dl,d,du,u);
379  // Solve solve state sensitivity system at current time step
380  linear_solve(ijv,dl,d,du,v);
381  }
382 
383  // Apply adjoint PDE Jacobian to a vector
384  void apply_adjoint_pde_jacobian(std::vector<Real> &ajv,
385  const std::vector<Real> &v,
386  const std::vector<Real> &u,
387  const std::vector<Real> &z) const {
388  // Fill jvp
389  for (int i = 0; i < nx_; i++) {
390  ajv[i] = nu_/dx_*2.0*v[i];
391  if ( i > 0 ) {
392  ajv[i] += -nu_/dx_*v[i-1]-nl_*(u[i-1]/6.0*v[i]
393  -(u[i-1]+2.0*u[i])/6.0*v[i-1]);
394  }
395  if ( i < nx_-1 ) {
396  ajv[i] += -nu_/dx_*v[i+1]+nl_*(u[i+1]/6.0*v[i]
397  -(u[i+1]+2.0*u[i])/6.0*v[i+1]);
398  }
399  }
400  ajv[ 0] -= nl_*u0_/6.0*v[0];
401  ajv[nx_-1] += nl_*u1_/6.0*v[nx_-1];
402  }
403 
404  // Apply inverse adjoint PDE Jacobian to a vector
405  void apply_inverse_adjoint_pde_jacobian(std::vector<Real> &iajv,
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> du(nx_-1,0.0);
412  std::vector<Real> dl(nx_-1,0.0);
413  compute_pde_jacobian(dl,d,du,u);
414  // Solve solve adjoint system at current time step
415  linear_solve(iajv,dl,d,du,v,true);
416  }
417 
418  /***************************************************************************/
419  /*********************** CONTROL JACOBIAN FUNCTIONS ************************/
420  /***************************************************************************/
421  // Apply control Jacobian to a vector
422  void apply_control_jacobian(std::vector<Real> &jv,
423  const std::vector<Real> &v,
424  const std::vector<Real> &u,
425  const std::vector<Real> &z) const {
426  for (int i=0; i<nx_; i++) {
427  // Contribution from control
428  jv[i] = -dx_/6.0*(v[i]+4.0*v[i+1]+v[i+2]);
429  }
430  }
431 
432  // Apply adjoint control Jacobian to a vector
433  void apply_adjoint_control_jacobian(std::vector<Real> &jv,
434  const std::vector<Real> &v,
435  const std::vector<Real> &u,
436  const std::vector<Real> &z) const {
437  for (int i=0; i<nx_+2; i++) {
438  if ( i == 0 ) {
439  jv[i] = -dx_/6.0*v[i];
440  }
441  else if ( i == 1 ) {
442  jv[i] = -dx_/6.0*(4.0*v[i-1]+v[i]);
443  }
444  else if ( i == nx_ ) {
445  jv[i] = -dx_/6.0*(4.0*v[i-1]+v[i-2]);
446  }
447  else if ( i == nx_+1 ) {
448  jv[i] = -dx_/6.0*v[i-2];
449  }
450  else {
451  jv[i] = -dx_/6.0*(v[i-2]+4.0*v[i-1]+v[i]);
452  }
453  }
454  }
455 
456  /***************************************************************************/
457  /*********************** AJDOINT HESSIANS **********************************/
458  /***************************************************************************/
459  void apply_adjoint_pde_hessian(std::vector<Real> &ahwv,
460  const std::vector<Real> &w,
461  const std::vector<Real> &v,
462  const std::vector<Real> &u,
463  const std::vector<Real> &z) const {
464  for (int i=0; i<nx_; i++) {
465  // Contribution from nonlinear term
466  ahwv[i] = 0.0;
467  if (i<nx_-1){
468  ahwv[i] += (w[i]*v[i+1] - w[i+1]*(2.0*v[i]+v[i+1]))/6.0;
469  }
470  if (i>0) {
471  ahwv[i] += (w[i-1]*(v[i-1]+2.0*v[i]) - w[i]*v[i-1])/6.0;
472  }
473  }
474  //ahwv.assign(u.size(),0.0);
475  }
476  void apply_adjoint_pde_control_hessian(std::vector<Real> &ahwv,
477  const std::vector<Real> &w,
478  const std::vector<Real> &v,
479  const std::vector<Real> &u,
480  const std::vector<Real> &z) {
481  ahwv.assign(u.size(),0.0);
482  }
483  void apply_adjoint_control_pde_hessian(std::vector<Real> &ahwv,
484  const std::vector<Real> &w,
485  const std::vector<Real> &v,
486  const std::vector<Real> &u,
487  const std::vector<Real> &z) {
488  ahwv.assign(z.size(),0.0);
489  }
490  void apply_adjoint_control_hessian(std::vector<Real> &ahwv,
491  const std::vector<Real> &w,
492  const std::vector<Real> &v,
493  const std::vector<Real> &u,
494  const std::vector<Real> &z) {
495  ahwv.assign(z.size(),0.0);
496  }
497 };
498 
499 template<class Real>
500 class L2VectorPrimal : public ROL::StdVector<Real> {
501 private:
502  ROL::Ptr<BurgersFEM<Real>> fem_;
503  mutable ROL::Ptr<L2VectorDual<Real>> dual_vec_;
504  mutable bool isDualInitialized_;
505 
506 public:
507  L2VectorPrimal(const ROL::Ptr<std::vector<Real>> & vec,
508  const ROL::Ptr<BurgersFEM<Real>> &fem)
509  : ROL::StdVector<Real>(vec),
510  fem_(fem), dual_vec_(ROL::nullPtr), isDualInitialized_(false) {}
511 
512  Real dot( const ROL::Vector<Real> &x ) const override {
513  const L2VectorPrimal & ex = dynamic_cast<const L2VectorPrimal&>(x);
514  const std::vector<Real>& xval = *ex.getVector();
515  const std::vector<Real>& yval = *ROL::StdVector<Real>::getVector();
516  return fem_->compute_L2_dot(xval,yval);
517  }
518 
519  ROL::Ptr<ROL::Vector<Real>> clone() const override {
520  size_t dimension = ROL::StdVector<Real>::getVector()->size();
521  return ROL::makePtr<L2VectorPrimal>(
522  ROL::makePtr<std::vector<Real>>(dimension,Real(0)),fem_);
523  }
524 
525  const ROL::Vector<Real>& dual() const override {
526  size_t dimension = ROL::StdVector<Real>::getVector()->size();
527  if ( !isDualInitialized_ ) {
528  dual_vec_ = ROL::makePtr<L2VectorDual<Real>>(
529  ROL::makePtr<std::vector<Real>>(dimension,Real(0)),fem_);
530  isDualInitialized_ = true;
531  }
532  fem_->apply_mass(*ROL::constPtrCast<std::vector<Real>>(dual_vec_->getVector()),
534  return *dual_vec_;
535  }
536 };
537 
538 template<class Real>
539 class L2VectorDual : public ROL::StdVector<Real> {
540 private:
541  ROL::Ptr<BurgersFEM<Real>> fem_;
542  mutable ROL::Ptr<L2VectorPrimal<Real>> dual_vec_;
543  mutable bool isDualInitialized_;
544 
545 public:
546  L2VectorDual(const ROL::Ptr<std::vector<Real>> & vec,
547  const ROL::Ptr<BurgersFEM<Real>> &fem)
548  : ROL::StdVector<Real>(vec),
549  fem_(fem), dual_vec_(ROL::nullPtr), isDualInitialized_(false) {}
550 
551  Real dot( const ROL::Vector<Real> &x ) const override {
552  const L2VectorDual & ex = dynamic_cast<const L2VectorDual&>(x);
553  const std::vector<Real>& xval = *ex.getVector();
554  const std::vector<Real>& yval = *ROL::StdVector<Real>::getVector();
555  size_t dimension = yval.size();
556  std::vector<Real> Mx(dimension,Real(0));
557  fem_->apply_inverse_mass(Mx,xval);
558  Real val(0);
559  for (size_t i = 0; i < dimension; i++) {
560  val += Mx[i]*yval[i];
561  }
562  return val;
563  }
564 
565  ROL::Ptr<ROL::Vector<Real>> clone() const override {
566  size_t dimension = ROL::StdVector<Real>::getVector()->size();
567  return ROL::makePtr<L2VectorDual>(
568  ROL::makePtr<std::vector<Real>>(dimension,Real(0)),fem_);
569  }
570 
571  const ROL::Vector<Real>& dual() const override {
572  size_t dimension = ROL::StdVector<Real>::getVector()->size();
573  if ( !isDualInitialized_ ) {
574  dual_vec_ = ROL::makePtr<L2VectorPrimal<Real>>(
575  ROL::makePtr<std::vector<Real>>(dimension,Real(0)),fem_);
576  isDualInitialized_ = true;
577  }
578  fem_->apply_inverse_mass(*ROL::constPtrCast<std::vector<Real>>(dual_vec_->getVector()),
580  return *dual_vec_;
581  }
582 };
583 
584 template<class Real>
585 class H1VectorPrimal : public ROL::StdVector<Real> {
586 private:
587  ROL::Ptr<BurgersFEM<Real>> fem_;
588  mutable ROL::Ptr<H1VectorDual<Real>> dual_vec_;
589  mutable bool isDualInitialized_;
590 
591 public:
592  H1VectorPrimal(const ROL::Ptr<std::vector<Real>> & vec,
593  const ROL::Ptr<BurgersFEM<Real>> &fem)
594  : ROL::StdVector<Real>(vec),
595  fem_(fem), dual_vec_(ROL::nullPtr), isDualInitialized_(false) {}
596 
597  Real dot( const ROL::Vector<Real> &x ) const override {
598  const H1VectorPrimal & ex = dynamic_cast<const H1VectorPrimal&>(x);
599  const std::vector<Real>& xval = *ex.getVector();
600  const std::vector<Real>& yval = *ROL::StdVector<Real>::getVector();
601  return fem_->compute_H1_dot(xval,yval);
602  }
603 
604  ROL::Ptr<ROL::Vector<Real>> clone() const override {
605  size_t dimension = ROL::StdVector<Real>::getVector()->size();
606  return ROL::makePtr<H1VectorPrimal>(
607  ROL::makePtr<std::vector<Real>>(dimension,Real(0)),fem_);
608  }
609 
610  const ROL::Vector<Real>& dual() const override {
611  size_t dimension = ROL::StdVector<Real>::getVector()->size();
612  if ( !isDualInitialized_ ) {
613  dual_vec_ = ROL::makePtr<H1VectorDual<Real>>(
614  ROL::makePtr<std::vector<Real>>(dimension,Real(0)),fem_);
615  isDualInitialized_ = true;
616  }
617  fem_->apply_H1(*ROL::constPtrCast<std::vector<Real>>(dual_vec_->getVector()),
619  return *dual_vec_;
620  }
621 };
622 
623 template<class Real>
624 class H1VectorDual : public ROL::StdVector<Real> {
625 private:
626  ROL::Ptr<BurgersFEM<Real>> fem_;
627  mutable ROL::Ptr<H1VectorPrimal<Real>> dual_vec_;
628  mutable bool isDualInitialized_;
629 
630 public:
631  H1VectorDual(const ROL::Ptr<std::vector<Real>> & vec,
632  const ROL::Ptr<BurgersFEM<Real>> &fem)
633  : ROL::StdVector<Real>(vec),
634  fem_(fem), dual_vec_(ROL::nullPtr), isDualInitialized_(false) {}
635 
636  Real dot( const ROL::Vector<Real> &x ) const override {
637  const H1VectorDual & ex = dynamic_cast<const H1VectorDual&>(x);
638  const std::vector<Real>& xval = *ex.getVector();
639  const std::vector<Real>& yval = *ROL::StdVector<Real>::getVector();
640  size_t dimension = yval.size();
641  std::vector<Real> Mx(dimension,0.0);
642  fem_->apply_inverse_H1(Mx,xval);
643  Real val(0);
644  for (size_t i = 0; i < dimension; i++) {
645  val += Mx[i]*yval[i];
646  }
647  return val;
648  }
649 
650  ROL::Ptr<ROL::Vector<Real>> clone() const override {
651  size_t dimension = ROL::StdVector<Real>::getVector()->size();
652  return ROL::makePtr<H1VectorDual>(
653  ROL::makePtr<std::vector<Real>>(dimension,Real(0)),fem_);
654  }
655 
656  const ROL::Vector<Real>& dual() const override {
657  size_t dimension = ROL::StdVector<Real>::getVector()->size();
658  if ( !isDualInitialized_ ) {
659  dual_vec_ = ROL::makePtr<H1VectorPrimal<Real>>(
660  ROL::makePtr<std::vector<Real>>(dimension,Real(0)),fem_);
661  isDualInitialized_ = true;
662  }
663  fem_->apply_inverse_H1(*ROL::constPtrCast<std::vector<Real>>(dual_vec_->getVector()),
665  return *dual_vec_;
666  }
667 };
668 
669 template<class Real>
671 private:
672 
675 
678 
681 
682  ROL::Ptr<BurgersFEM<Real> > fem_;
684 
685 public:
687  const bool useHessian = true)
688  : fem_(fem), useHessian_(useHessian) {}
689 
691  const ROL::Vector<Real> &z, Real &tol) {
692  ROL::Ptr<std::vector<Real> > cp =
693  dynamic_cast<PrimalConstraintVector&>(c).getVector();
694  ROL::Ptr<const std::vector<Real> > up =
695  dynamic_cast<const PrimalStateVector&>(u).getVector();
696  ROL::Ptr<const std::vector<Real> > zp =
697  dynamic_cast<const PrimalControlVector&>(z).getVector();
698 
699  fem_->compute_residual(*cp,*up,*zp);
700  }
701 
703  const ROL::Vector<Real> &z, Real &tol) {
704  ROL::Ptr<std::vector<Real> > jvp =
705  dynamic_cast<PrimalConstraintVector&>(jv).getVector();
706  ROL::Ptr<const std::vector<Real> > vp =
707  dynamic_cast<const PrimalStateVector&>(v).getVector();
708  ROL::Ptr<const std::vector<Real> > up =
709  dynamic_cast<const PrimalStateVector&>(u).getVector();
710  ROL::Ptr<const std::vector<Real> > zp =
711  dynamic_cast<const PrimalControlVector&>(z).getVector();
712 
713  fem_->apply_pde_jacobian(*jvp,*vp,*up,*zp);
714  }
715 
717  const ROL::Vector<Real> &z, Real &tol) {
718  ROL::Ptr<std::vector<Real> > jvp =
719  dynamic_cast<PrimalConstraintVector&>(jv).getVector();
720  ROL::Ptr<const std::vector<Real> > vp =
721  dynamic_cast<const PrimalControlVector&>(v).getVector();
722  ROL::Ptr<const std::vector<Real> > up =
723  dynamic_cast<const PrimalStateVector&>(u).getVector();
724  ROL::Ptr<const std::vector<Real> > zp =
725  dynamic_cast<const PrimalControlVector&>(z).getVector();
726 
727  fem_->apply_control_jacobian(*jvp,*vp,*up,*zp);
728  }
729 
731  const ROL::Vector<Real> &z, Real &tol) {
732  ROL::Ptr<std::vector<Real> > ijvp =
733  dynamic_cast<PrimalStateVector&>(ijv).getVector();
734  ROL::Ptr<const std::vector<Real> > vp =
735  dynamic_cast<const PrimalConstraintVector&>(v).getVector();
736  ROL::Ptr<const std::vector<Real> > up =
737  dynamic_cast<const PrimalStateVector&>(u).getVector();
738  ROL::Ptr<const std::vector<Real> > zp =
739  dynamic_cast<const PrimalControlVector&>(z).getVector();
740 
741  fem_->apply_inverse_pde_jacobian(*ijvp,*vp,*up,*zp);
742  }
743 
745  const ROL::Vector<Real> &z, Real &tol) {
746  ROL::Ptr<std::vector<Real> > jvp =
747  dynamic_cast<DualStateVector&>(ajv).getVector();
748  ROL::Ptr<const std::vector<Real> > vp =
749  dynamic_cast<const DualConstraintVector&>(v).getVector();
750  ROL::Ptr<const std::vector<Real> > up =
751  dynamic_cast<const PrimalStateVector&>(u).getVector();
752  ROL::Ptr<const std::vector<Real> > zp =
753  dynamic_cast<const PrimalControlVector&>(z).getVector();
754 
755  fem_->apply_adjoint_pde_jacobian(*jvp,*vp,*up,*zp);
756  }
757 
759  const ROL::Vector<Real> &z, Real &tol) {
760  ROL::Ptr<std::vector<Real> > jvp =
761  dynamic_cast<DualControlVector&>(jv).getVector();
762  ROL::Ptr<const std::vector<Real> > vp =
763  dynamic_cast<const DualConstraintVector&>(v).getVector();
764  ROL::Ptr<const std::vector<Real> > up =
765  dynamic_cast<const PrimalStateVector&>(u).getVector();
766  ROL::Ptr<const std::vector<Real> > zp =
767  dynamic_cast<const PrimalControlVector&>(z).getVector();
768 
769  fem_->apply_adjoint_control_jacobian(*jvp,*vp,*up,*zp);
770  }
771 
773  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
774  ROL::Ptr<std::vector<Real> > iajvp =
775  dynamic_cast<DualConstraintVector&>(iajv).getVector();
776  ROL::Ptr<const std::vector<Real> > vp =
777  dynamic_cast<const DualStateVector&>(v).getVector();
778  ROL::Ptr<const std::vector<Real> > up =
779  dynamic_cast<const PrimalStateVector&>(u).getVector();
780  ROL::Ptr<const std::vector<Real> > zp =
781  dynamic_cast<const PrimalControlVector&>(z).getVector();
782 
783  fem_->apply_inverse_adjoint_pde_jacobian(*iajvp,*vp,*up,*zp);
784  }
785 
787  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
788  if ( useHessian_ ) {
789  ROL::Ptr<std::vector<Real> > ahwvp =
790  dynamic_cast<DualStateVector&>(ahwv).getVector();
791  ROL::Ptr<const std::vector<Real> > wp =
792  dynamic_cast<const DualConstraintVector&>(w).getVector();
793  ROL::Ptr<const std::vector<Real> > vp =
794  dynamic_cast<const PrimalStateVector&>(v).getVector();
795  ROL::Ptr<const std::vector<Real> > up =
796  dynamic_cast<const PrimalStateVector&>(u).getVector();
797  ROL::Ptr<const std::vector<Real> > zp =
798  dynamic_cast<const PrimalControlVector&>(z).getVector();
799 
800  fem_->apply_adjoint_pde_hessian(*ahwvp,*wp,*vp,*up,*zp);
801  }
802  else {
803  ahwv.zero();
804  }
805  }
806 
808  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
809  if ( useHessian_ ) {
810  ROL::Ptr<std::vector<Real> > ahwvp =
811  dynamic_cast<DualControlVector&>(ahwv).getVector();
812  ROL::Ptr<const std::vector<Real> > wp =
813  dynamic_cast<const DualConstraintVector&>(w).getVector();
814  ROL::Ptr<const std::vector<Real> > vp =
815  dynamic_cast<const PrimalStateVector&>(v).getVector();
816  ROL::Ptr<const std::vector<Real> > up =
817  dynamic_cast<const PrimalStateVector&>(u).getVector();
818  ROL::Ptr<const std::vector<Real> > zp =
819  dynamic_cast<const PrimalControlVector&>(z).getVector();
820 
821  fem_->apply_adjoint_control_pde_hessian(*ahwvp,*wp,*vp,*up,*zp);
822  }
823  else {
824  ahwv.zero();
825  }
826  }
828  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
829  if ( useHessian_ ) {
830  ROL::Ptr<std::vector<Real> > ahwvp =
831  dynamic_cast<DualStateVector&>(ahwv).getVector();
832  ROL::Ptr<const std::vector<Real> > wp =
833  dynamic_cast<const DualConstraintVector&>(w).getVector();
834  ROL::Ptr<const std::vector<Real> > vp =
835  dynamic_cast<const PrimalControlVector&>(v).getVector();
836  ROL::Ptr<const std::vector<Real> > up =
837  dynamic_cast<const PrimalStateVector&>(u).getVector();
838  ROL::Ptr<const std::vector<Real> > zp =
839  dynamic_cast<const PrimalControlVector&>(z).getVector();
840 
841  fem_->apply_adjoint_pde_control_hessian(*ahwvp,*wp,*vp,*up,*zp);
842  }
843  else {
844  ahwv.zero();
845  }
846  }
848  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
849  if ( useHessian_ ) {
850  ROL::Ptr<std::vector<Real> > ahwvp =
851  dynamic_cast<DualControlVector&>(ahwv).getVector();
852  ROL::Ptr<const std::vector<Real> > wp =
853  dynamic_cast<const DualConstraintVector&>(w).getVector();
854  ROL::Ptr<const std::vector<Real> > vp =
855  dynamic_cast<const PrimalControlVector&>(v).getVector();
856  ROL::Ptr<const std::vector<Real> > up =
857  dynamic_cast<const PrimalStateVector&>(u).getVector();
858  ROL::Ptr<const std::vector<Real> > zp =
859  dynamic_cast<const PrimalControlVector&>(z).getVector();
860 
861  fem_->apply_adjoint_control_hessian(*ahwvp,*wp,*vp,*up,*zp);
862  }
863  else {
864  ahwv.zero();
865  }
866  }
867 };
H1VectorPrimal< Real > DualConstraintVector
Definition: test_04.hpp:680
BurgersFEM(int nx=128, Real nl=1.0, Real cH1=1.0, Real cL2=1.0)
Definition: test_04.hpp:97
Real dx_
Definition: test_04.hpp:38
Real compute_H1_dot(const std::vector< Real > &x, const std::vector< Real > &y) const
Definition: test_04.hpp:214
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:744
Real cL2_
Definition: test_04.hpp:45
void apply_mass(std::vector< Real > &Mu, const std::vector< Real > &u) const
Definition: test_04.hpp:147
void axpy(std::vector< Real > &out, const Real a, const std::vector< Real > &x, const std::vector< Real > &y) const
Definition: test_04.hpp:54
ROL::Ptr< ROL::Vector< Real > > clone() const override
Clone to make a new (uninitialized) vector.
Definition: test_04.hpp:604
Real dot(const ROL::Vector< Real > &x) const override
Compute where .
Definition: test_04.hpp:636
bool isDualInitialized_
Definition: test_04.hpp:628
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:476
const ROL::Vector< Real > & dual() const override
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
Definition: test_04.hpp:525
int dimension() const
Return dimension of the vector space.
Definition: example_04.hpp:558
ROL::Ptr< L2VectorDual< Real > > dual_vec_
Definition: test_04.hpp:503
StdVector(const Ptr< std::vector< Real >> &std_vec)
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:807
Constraint_BurgersControl(const ROL::Ptr< BurgersFEM< Real > > &fem, const bool useHessian=true)
Definition: test_04.hpp:686
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:459
ROL::Ptr< BurgersFEM< Real > > fem_
Definition: test_04.hpp:626
bool isDualInitialized_
Definition: test_04.hpp:543
Real u0_
Definition: test_04.hpp:41
Contains definitions of custom data types in ROL.
Ptr< const std::vector< Element > > getVector() const
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:786
ROL::Ptr< const std::vector< Real > > getVector() const
Definition: example_04.hpp:543
ROL::Ptr< BurgersFEM< Real > > fem_
Definition: test_04.hpp:587
Real compute_H1_norm(const std::vector< Real > &r) const
Definition: test_04.hpp:234
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:772
const ROL::Vector< Real > & dual() const override
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
Definition: test_04.hpp:656
ROL::Ptr< ROL::Vector< Real > > clone() const override
Clone to make a new (uninitialized) vector.
Definition: test_04.hpp:519
bool isDualInitialized_
Definition: test_04.hpp:504
virtual void zero()
Set to zero vector.
Definition: ROL_Vector.hpp:133
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:46
H1VectorDual< Real > PrimalConstraintVector
Definition: test_04.hpp:679
L2VectorPrimal(const ROL::Ptr< std::vector< Real >> &vec, const ROL::Ptr< BurgersFEM< Real >> &fem)
Definition: test_04.hpp:507
Real dot(const ROL::Vector< Real > &x) const override
Compute where .
Definition: test_04.hpp:597
int num_dof(void) const
Definition: test_04.hpp:112
ROL::Ptr< ROL::Vector< Real > > clone() const override
Clone to make a new (uninitialized) vector.
Definition: test_04.hpp:565
H1VectorPrimal< Real > PrimalStateVector
Definition: test_04.hpp:673
Real mesh_spacing(void) const
Definition: test_04.hpp:116
int dimension() const
Return dimension of the vector space.
Definition: example_04.hpp:737
void test_inverse_mass(std::ostream &outStream=std::cout)
Definition: test_04.hpp:177
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:351
H1VectorDual(const ROL::Ptr< std::vector< Real >> &vec, const ROL::Ptr< BurgersFEM< Real >> &fem)
Definition: test_04.hpp:631
ROL::Ptr< BurgersFEM< Real > > fem_
Definition: test_04.hpp:682
Real nl_
Definition: test_04.hpp:40
Real cH1_
Definition: test_04.hpp:44
void applyAdjointHessian_21(ROL::Vector< Real > &ahwv, const ROL::Vector< Real > &w, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Apply the simulation-space derivative of the adjoint of the constraint optimization-space Jacobian at...
Definition: test_04.hpp:827
Real u1_
Definition: test_04.hpp:42
Provides the ROL::Vector interface for scalar values, to be used, for example, with scalar constraint...
ROL::Ptr< const std::vector< Real > > getVector() const
Definition: example_04.hpp:722
ROL::Ptr< const std::vector< Real > > getVector() const
Definition: example_04.hpp:636
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:323
H1VectorDual< Real > DualStateVector
Definition: test_04.hpp:674
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:702
Real compute_L2_dot(const std::vector< Real > &x, const std::vector< Real > &y) const
Definition: test_04.hpp:124
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:66
L2VectorPrimal< Real > PrimalControlVector
Definition: test_04.hpp:676
void apply_inverse_mass(std::vector< Real > &Mu, const std::vector< Real > &u) const
Definition: test_04.hpp:164
L2VectorDual(const ROL::Ptr< std::vector< Real >> &vec, const ROL::Ptr< BurgersFEM< Real >> &fem)
Definition: test_04.hpp:546
ROL::Ptr< BurgersFEM< Real > > fem_
Definition: test_04.hpp:502
ROL::Ptr< const std::vector< Real > > getVector() const
Definition: example_04.hpp:815
ROL::Ptr< BurgersFEM< Real > > fem_
Definition: test_04.hpp:541
void update(std::vector< Real > &u, const std::vector< Real > &s, const Real alpha=1.0) const
Definition: test_04.hpp:48
const ROL::Vector< Real > & dual() const override
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
Definition: test_04.hpp:571
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:490
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:730
void compute_residual(std::vector< Real > &r, const std::vector< Real > &u, const std::vector< Real > &z) const
Definition: test_04.hpp:287
Real nu_
Definition: test_04.hpp:39
void set_problem_data(const Real nu, const Real f, const Real u0, const Real u1)
Definition: test_04.hpp:105
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:716
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:422
Real dot(const ROL::Vector< Real > &x) const override
Compute where .
Definition: test_04.hpp:551
void test_inverse_H1(std::ostream &outStream=std::cout)
Definition: test_04.hpp:266
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:433
H1VectorPrimal(const ROL::Ptr< std::vector< Real >> &vec, const ROL::Ptr< BurgersFEM< Real >> &fem)
Definition: test_04.hpp:592
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:370
L2VectorDual< Real > DualControlVector
Definition: test_04.hpp:677
Real f_
Definition: test_04.hpp:43
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:847
void apply_H1(std::vector< Real > &Mu, const std::vector< Real > &u) const
Definition: test_04.hpp:239
void apply_inverse_H1(std::vector< Real > &Mu, const std::vector< Real > &u) const
Definition: test_04.hpp:258
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:384
ROL::Ptr< ROL::Vector< Real > > clone() const override
Clone to make a new (uninitialized) vector.
Definition: test_04.hpp:650
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:405
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:483
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:758
int dimension() const
Return dimension of the vector space.
Definition: example_04.hpp:830
Defines the constraint operator interface for simulation-based optimization.
ROL::Ptr< H1VectorPrimal< Real > > dual_vec_
Definition: test_04.hpp:627
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:690
void scale(std::vector< Real > &u, const Real alpha=0.0) const
Definition: test_04.hpp:60
bool isDualInitialized_
Definition: test_04.hpp:589
constexpr auto dim
const ROL::Vector< Real > & dual() const override
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
Definition: test_04.hpp:610
ROL::Ptr< H1VectorDual< Real > > dual_vec_
Definition: test_04.hpp:588
Real dot(const ROL::Vector< Real > &x) const override
Compute where .
Definition: test_04.hpp:512
int dimension() const
Return dimension of the vector space.
Definition: example_04.hpp:651
Real compute_L2_norm(const std::vector< Real > &r) const
Definition: test_04.hpp:142
ROL::Ptr< L2VectorPrimal< Real > > dual_vec_
Definition: test_04.hpp:542