ROL
example_04.hpp
Go to the documentation of this file.
1// @HEADER
2// ************************************************************************
3//
4// Rapid Optimization Library (ROL) Package
5// Copyright (2014) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact lead developers:
38// Drew Kouri (dpkouri@sandia.gov) and
39// Denis Ridzal (dridzal@sandia.gov)
40//
41// ************************************************************************
42// @HEADER
43
49#include "ROL_Types.hpp"
50#include "ROL_Vector.hpp"
54
55#include "Teuchos_LAPACK.hpp"
56
57template<class Real>
58class L2VectorPrimal;
59
60template<class Real>
61class L2VectorDual;
62
63template<class Real>
64class H1VectorPrimal;
65
66template<class Real>
67class H1VectorDual;
68
69template<class Real>
70class BurgersFEM {
71private:
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
82private:
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
131public:
132 BurgersFEM(int nx = 128, Real nu = 1.e-2, Real nl = 1.0,
133 Real u0 = 1.0, Real u1 = 0.0, Real f = 0.0,
134 Real cH1 = 1.0, Real cL2 = 1.0)
135 : nx_(nx), dx_(1.0/((Real)nx+1.0)),
136 nu_(nu), nl_(nl), u0_(u0), u1_(u1), f_(f),
137 cH1_(cH1), cL2_(cL2) {}
138
139 int num_dof(void) const {
140 return nx_;
141 }
142
143 Real mesh_spacing(void) const {
144 return dx_;
145 }
146
147 /***************************************************************************/
148 /*********************** L2 INNER PRODUCT FUNCTIONS ************************/
149 /***************************************************************************/
150 // Compute L2 inner product
151 Real compute_L2_dot(const std::vector<Real> &x, const std::vector<Real> &y) const {
152 Real ip = 0.0;
153 Real c = (((int)x.size()==nx_) ? 4.0 : 2.0);
154 for (unsigned i=0; i<x.size(); i++) {
155 if ( i == 0 ) {
156 ip += dx_/6.0*(c*x[i] + x[i+1])*y[i];
157 }
158 else if ( i == x.size()-1 ) {
159 ip += dx_/6.0*(x[i-1] + c*x[i])*y[i];
160 }
161 else {
162 ip += dx_/6.0*(x[i-1] + 4.0*x[i] + x[i+1])*y[i];
163 }
164 }
165 return ip;
166 }
167
168 // compute L2 norm
169 Real compute_L2_norm(const std::vector<Real> &r) const {
170 return std::sqrt(compute_L2_dot(r,r));
171 }
172
173 // Apply L2 Reisz operator
174 void apply_mass(std::vector<Real> &Mu, const std::vector<Real> &u ) const {
175 Mu.resize(u.size(),0.0);
176 Real c = (((int)u.size()==nx_) ? 4.0 : 2.0);
177 for (unsigned i=0; i<u.size(); i++) {
178 if ( i == 0 ) {
179 Mu[i] = dx_/6.0*(c*u[i] + u[i+1]);
180 }
181 else if ( i == u.size()-1 ) {
182 Mu[i] = dx_/6.0*(u[i-1] + c*u[i]);
183 }
184 else {
185 Mu[i] = dx_/6.0*(u[i-1] + 4.0*u[i] + u[i+1]);
186 }
187 }
188 }
189
190 // Apply L2 inverse Reisz operator
191 void apply_inverse_mass(std::vector<Real> &Mu, const std::vector<Real> &u) const {
192 unsigned nx = u.size();
193 // Build mass matrix
194 std::vector<Real> dl(nx-1,dx_/6.0);
195 std::vector<Real> d(nx,2.0*dx_/3.0);
196 std::vector<Real> du(nx-1,dx_/6.0);
197 if ( (int)nx != nx_ ) {
198 d[ 0] = dx_/3.0;
199 d[nx-1] = dx_/3.0;
200 }
201 linear_solve(Mu,dl,d,du,u);
202 }
203
204 void test_inverse_mass(std::ostream &outStream = std::cout) {
205 // State Mass Matrix
206 std::vector<Real> u(nx_,0.0), Mu(nx_,0.0), iMMu(nx_,0.0), diff(nx_,0.0);
207 for (int i = 0; i < nx_; i++) {
208 u[i] = 2.0*(Real)rand()/(Real)RAND_MAX - 1.0;
209 }
210 apply_mass(Mu,u);
211 apply_inverse_mass(iMMu,Mu);
212 axpy(diff,-1.0,iMMu,u);
213 Real error = compute_L2_norm(diff);
214 Real normu = compute_L2_norm(u);
215 outStream << "Test Inverse State Mass Matrix\n";
216 outStream << " ||u - inv(M)Mu|| = " << error << "\n";
217 outStream << " ||u|| = " << normu << "\n";
218 outStream << " Relative Error = " << error/normu << "\n";
219 outStream << "\n";
220 // Control Mass Matrix
221 u.resize(nx_+2,0.0); Mu.resize(nx_+2,0.0); iMMu.resize(nx_+2,0.0); diff.resize(nx_+2,0.0);
222 for (int i = 0; i < nx_+2; i++) {
223 u[i] = 2.0*(Real)rand()/(Real)RAND_MAX - 1.0;
224 }
225 apply_mass(Mu,u);
226 apply_inverse_mass(iMMu,Mu);
227 axpy(diff,-1.0,iMMu,u);
228 error = compute_L2_norm(diff);
229 normu = compute_L2_norm(u);
230 outStream << "Test Inverse Control Mass Matrix\n";
231 outStream << " ||z - inv(M)Mz|| = " << error << "\n";
232 outStream << " ||z|| = " << normu << "\n";
233 outStream << " Relative Error = " << error/normu << "\n";
234 outStream << "\n";
235 }
236
237 /***************************************************************************/
238 /*********************** H1 INNER PRODUCT FUNCTIONS ************************/
239 /***************************************************************************/
240 // Compute H1 inner product
241 Real compute_H1_dot(const std::vector<Real> &x, const std::vector<Real> &y) const {
242 Real ip = 0.0;
243 for (int i=0; i<nx_; i++) {
244 if ( i == 0 ) {
245 ip += cL2_*dx_/6.0*(4.0*x[i] + x[i+1])*y[i]; // Mass term
246 ip += cH1_*(2.0*x[i] - x[i+1])/dx_*y[i]; // Stiffness term
247 }
248 else if ( i == nx_-1 ) {
249 ip += cL2_*dx_/6.0*(x[i-1] + 4.0*x[i])*y[i]; // Mass term
250 ip += cH1_*(2.0*x[i] - x[i-1])/dx_*y[i]; // Stiffness term
251 }
252 else {
253 ip += cL2_*dx_/6.0*(x[i-1] + 4.0*x[i] + x[i+1])*y[i]; // Mass term
254 ip += cH1_*(2.0*x[i] - x[i-1] - x[i+1])/dx_*y[i]; // Stiffness term
255 }
256 }
257 return ip;
258 }
259
260 // compute H1 norm
261 Real compute_H1_norm(const std::vector<Real> &r) const {
262 return std::sqrt(compute_H1_dot(r,r));
263 }
264
265 // Apply H2 Reisz operator
266 void apply_H1(std::vector<Real> &Mu, const std::vector<Real> &u ) const {
267 Mu.resize(nx_,0.0);
268 for (int i=0; i<nx_; i++) {
269 if ( i == 0 ) {
270 Mu[i] = cL2_*dx_/6.0*(4.0*u[i] + u[i+1])
271 + cH1_*(2.0*u[i] - u[i+1])/dx_;
272 }
273 else if ( i == nx_-1 ) {
274 Mu[i] = cL2_*dx_/6.0*(u[i-1] + 4.0*u[i])
275 + cH1_*(2.0*u[i] - u[i-1])/dx_;
276 }
277 else {
278 Mu[i] = cL2_*dx_/6.0*(u[i-1] + 4.0*u[i] + u[i+1])
279 + cH1_*(2.0*u[i] - u[i-1] - u[i+1])/dx_;
280 }
281 }
282 }
283
284 // Apply H1 inverse Reisz operator
285 void apply_inverse_H1(std::vector<Real> &Mu, const std::vector<Real> &u) const {
286 // Build mass matrix
287 std::vector<Real> dl(nx_-1,cL2_*dx_/6.0 - cH1_/dx_);
288 std::vector<Real> d(nx_,2.0*(cL2_*dx_/3.0 + cH1_/dx_));
289 std::vector<Real> du(nx_-1,cL2_*dx_/6.0 - cH1_/dx_);
290 linear_solve(Mu,dl,d,du,u);
291 }
292
293 void test_inverse_H1(std::ostream &outStream = std::cout) {
294 std::vector<Real> u(nx_,0.0), Mu(nx_,0.0), iMMu(nx_,0.0), diff(nx_,0.0);
295 for (int i = 0; i < nx_; i++) {
296 u[i] = 2.0*(Real)rand()/(Real)RAND_MAX - 1.0;
297 }
298 apply_H1(Mu,u);
299 apply_inverse_H1(iMMu,Mu);
300 axpy(diff,-1.0,iMMu,u);
301 Real error = compute_H1_norm(diff);
302 Real normu = compute_H1_norm(u);
303 outStream << "Test Inverse State H1 Matrix\n";
304 outStream << " ||u - inv(M)Mu|| = " << error << "\n";
305 outStream << " ||u|| = " << normu << "\n";
306 outStream << " Relative Error = " << error/normu << "\n";
307 outStream << "\n";
308 }
309
310 /***************************************************************************/
311 /*********************** PDE RESIDUAL AND SOLVE ****************************/
312 /***************************************************************************/
313 // Compute current PDE residual
314 void compute_residual(std::vector<Real> &r, const std::vector<Real> &u,
315 const std::vector<Real> &z) const {
316 r.clear();
317 r.resize(nx_,0.0);
318 for (int i=0; i<nx_; i++) {
319 // Contribution from stiffness term
320 if ( i==0 ) {
321 r[i] = nu_/dx_*(2.0*u[i]-u[i+1]);
322 }
323 else if (i==nx_-1) {
324 r[i] = nu_/dx_*(2.0*u[i]-u[i-1]);
325 }
326 else {
327 r[i] = nu_/dx_*(2.0*u[i]-u[i-1]-u[i+1]);
328 }
329 // Contribution from nonlinear term
330 if (i<nx_-1){
331 r[i] += nl_*u[i+1]*(u[i]+u[i+1])/6.0;
332 }
333 if (i>0) {
334 r[i] -= nl_*u[i-1]*(u[i-1]+u[i])/6.0;
335 }
336 // Contribution from control
337 r[i] -= dx_/6.0*(z[i]+4.0*z[i+1]+z[i+2]);
338 // Contribution from right-hand side
339 r[i] -= dx_*f_;
340 }
341 // Contribution from Dirichlet boundary terms
342 r[0] -= nl_*(u0_*u[ 0]/6.0 + u0_*u0_/6.0) + nu_*u0_/dx_;
343 r[nx_-1] += nl_*(u1_*u[nx_-1]/6.0 + u1_*u1_/6.0) - nu_*u1_/dx_;
344 }
345
346 /***************************************************************************/
347 /*********************** PDE JACOBIAN FUNCTIONS ****************************/
348 /***************************************************************************/
349 // Build PDE Jacobian trigiagonal matrix
350 void compute_pde_jacobian(std::vector<Real> &dl, // Lower diagonal
351 std::vector<Real> &d, // Diagonal
352 std::vector<Real> &du, // Upper diagonal
353 const std::vector<Real> &u) const { // State variable
354 // Get Diagonal and Off-Diagonal Entries of linear PDE Jacobian
355 d.clear();
356 d.resize(nx_,nu_*2.0/dx_);
357 dl.clear();
358 dl.resize(nx_-1,-nu_/dx_);
359 du.clear();
360 du.resize(nx_-1,-nu_/dx_);
361 // Contribution from nonlinearity
362 for (int i=0; i<nx_; i++) {
363 if (i<nx_-1) {
364 dl[i] += nl_*(-2.0*u[i]-u[i+1])/6.0;
365 d[i] += nl_*u[i+1]/6.0;
366 }
367 if (i>0) {
368 d[i] -= nl_*u[i-1]/6.0;
369 du[i-1] += nl_*(u[i-1]+2.0*u[i])/6.0;
370 }
371 }
372 // Contribution from Dirichlet boundary conditions
373 d[ 0] -= nl_*u0_/6.0;
374 d[nx_-1] += nl_*u1_/6.0;
375 }
376
377 // Apply PDE Jacobian to a vector
378 void apply_pde_jacobian(std::vector<Real> &jv,
379 const std::vector<Real> &v,
380 const std::vector<Real> &u,
381 const std::vector<Real> &z) const {
382 // Fill jv
383 for (int i = 0; i < nx_; i++) {
384 jv[i] = nu_/dx_*2.0*v[i];
385 if ( i > 0 ) {
386 jv[i] += -nu_/dx_*v[i-1]-nl_*(u[i-1]/6.0*v[i]+(u[i]+2.0*u[i-1])/6.0*v[i-1]);
387 }
388 if ( i < nx_-1 ) {
389 jv[i] += -nu_/dx_*v[i+1]+nl_*(u[i+1]/6.0*v[i]+(u[i]+2.0*u[i+1])/6.0*v[i+1]);
390 }
391 }
392 jv[ 0] -= nl_*u0_/6.0*v[0];
393 jv[nx_-1] += nl_*u1_/6.0*v[nx_-1];
394 }
395
396 // Apply inverse PDE Jacobian to a vector
397 void apply_inverse_pde_jacobian(std::vector<Real> &ijv,
398 const std::vector<Real> &v,
399 const std::vector<Real> &u,
400 const std::vector<Real> &z) const {
401 // Get PDE Jacobian
402 std::vector<Real> d(nx_,0.0);
403 std::vector<Real> dl(nx_-1,0.0);
404 std::vector<Real> du(nx_-1,0.0);
405 compute_pde_jacobian(dl,d,du,u);
406 // Solve solve state sensitivity system at current time step
407 linear_solve(ijv,dl,d,du,v);
408 }
409
410 // Apply adjoint PDE Jacobian to a vector
411 void apply_adjoint_pde_jacobian(std::vector<Real> &ajv,
412 const std::vector<Real> &v,
413 const std::vector<Real> &u,
414 const std::vector<Real> &z) const {
415 // Fill jvp
416 for (int i = 0; i < nx_; i++) {
417 ajv[i] = nu_/dx_*2.0*v[i];
418 if ( i > 0 ) {
419 ajv[i] += -nu_/dx_*v[i-1]-nl_*(u[i-1]/6.0*v[i]
420 -(u[i-1]+2.0*u[i])/6.0*v[i-1]);
421 }
422 if ( i < nx_-1 ) {
423 ajv[i] += -nu_/dx_*v[i+1]+nl_*(u[i+1]/6.0*v[i]
424 -(u[i+1]+2.0*u[i])/6.0*v[i+1]);
425 }
426 }
427 ajv[ 0] -= nl_*u0_/6.0*v[0];
428 ajv[nx_-1] += nl_*u1_/6.0*v[nx_-1];
429 }
430
431 // Apply inverse adjoint PDE Jacobian to a vector
432 void apply_inverse_adjoint_pde_jacobian(std::vector<Real> &iajv,
433 const std::vector<Real> &v,
434 const std::vector<Real> &u,
435 const std::vector<Real> &z) const {
436 // Get PDE Jacobian
437 std::vector<Real> d(nx_,0.0);
438 std::vector<Real> du(nx_-1,0.0);
439 std::vector<Real> dl(nx_-1,0.0);
440 compute_pde_jacobian(dl,d,du,u);
441 // Solve solve adjoint system at current time step
442 linear_solve(iajv,dl,d,du,v,true);
443 }
444
445 /***************************************************************************/
446 /*********************** CONTROL JACOBIAN FUNCTIONS ************************/
447 /***************************************************************************/
448 // Apply control Jacobian to a vector
449 void apply_control_jacobian(std::vector<Real> &jv,
450 const std::vector<Real> &v,
451 const std::vector<Real> &u,
452 const std::vector<Real> &z) const {
453 for (int i=0; i<nx_; i++) {
454 // Contribution from control
455 jv[i] = -dx_/6.0*(v[i]+4.0*v[i+1]+v[i+2]);
456 }
457 }
458
459 // Apply adjoint control Jacobian to a vector
460 void apply_adjoint_control_jacobian(std::vector<Real> &jv,
461 const std::vector<Real> &v,
462 const std::vector<Real> &u,
463 const std::vector<Real> &z) const {
464 for (int i=0; i<nx_+2; i++) {
465 if ( i == 0 ) {
466 jv[i] = -dx_/6.0*v[i];
467 }
468 else if ( i == 1 ) {
469 jv[i] = -dx_/6.0*(4.0*v[i-1]+v[i]);
470 }
471 else if ( i == nx_ ) {
472 jv[i] = -dx_/6.0*(4.0*v[i-1]+v[i-2]);
473 }
474 else if ( i == nx_+1 ) {
475 jv[i] = -dx_/6.0*v[i-2];
476 }
477 else {
478 jv[i] = -dx_/6.0*(v[i-2]+4.0*v[i-1]+v[i]);
479 }
480 }
481 }
482
483 /***************************************************************************/
484 /*********************** AJDOINT HESSIANS **********************************/
485 /***************************************************************************/
486 void apply_adjoint_pde_hessian(std::vector<Real> &ahwv,
487 const std::vector<Real> &w,
488 const std::vector<Real> &v,
489 const std::vector<Real> &u,
490 const std::vector<Real> &z) const {
491 for (int i=0; i<nx_; i++) {
492 // Contribution from nonlinear term
493 ahwv[i] = 0.0;
494 if (i<nx_-1){
495 ahwv[i] += (w[i]*v[i+1] - w[i+1]*(2.0*v[i]+v[i+1]))/6.0;
496 }
497 if (i>0) {
498 ahwv[i] += (w[i-1]*(v[i-1]+2.0*v[i]) - w[i]*v[i-1])/6.0;
499 }
500 }
501 //ahwv.assign(u.size(),0.0);
502 }
503 void apply_adjoint_pde_control_hessian(std::vector<Real> &ahwv,
504 const std::vector<Real> &w,
505 const std::vector<Real> &v,
506 const std::vector<Real> &u,
507 const std::vector<Real> &z) {
508 ahwv.assign(u.size(),0.0);
509 }
510 void apply_adjoint_control_pde_hessian(std::vector<Real> &ahwv,
511 const std::vector<Real> &w,
512 const std::vector<Real> &v,
513 const std::vector<Real> &u,
514 const std::vector<Real> &z) {
515 ahwv.assign(z.size(),0.0);
516 }
517 void apply_adjoint_control_hessian(std::vector<Real> &ahwv,
518 const std::vector<Real> &w,
519 const std::vector<Real> &v,
520 const std::vector<Real> &u,
521 const std::vector<Real> &z) {
522 ahwv.assign(z.size(),0.0);
523 }
524};
525
526template<class Real>
527class L2VectorPrimal : public ROL::Vector<Real> {
528private:
529 ROL::Ptr<std::vector<Real> > vec_;
530 ROL::Ptr<BurgersFEM<Real> > fem_;
531
532 mutable ROL::Ptr<L2VectorDual<Real> > dual_vec_;
533
534public:
535 L2VectorPrimal(const ROL::Ptr<std::vector<Real> > & vec,
536 const ROL::Ptr<BurgersFEM<Real> > &fem)
537 : vec_(vec), fem_(fem), dual_vec_(ROL::nullPtr) {}
538
539 void set( const ROL::Vector<Real> &x ) {
540 const L2VectorPrimal &ex = dynamic_cast<const L2VectorPrimal&>(x);
541 const std::vector<Real>& xval = *ex.getVector();
542 std::copy(xval.begin(),xval.end(),vec_->begin());
543 }
544
545 void plus( const ROL::Vector<Real> &x ) {
546 const L2VectorPrimal &ex = dynamic_cast<const L2VectorPrimal&>(x);
547 const std::vector<Real>& xval = *ex.getVector();
548 unsigned dimension = vec_->size();
549 for (unsigned i=0; i<dimension; i++) {
550 (*vec_)[i] += xval[i];
551 }
552 }
553
554 void scale( const Real alpha ) {
555 unsigned dimension = vec_->size();
556 for (unsigned i=0; i<dimension; i++) {
557 (*vec_)[i] *= alpha;
558 }
559 }
560
561 Real dot( const ROL::Vector<Real> &x ) const {
562 const L2VectorPrimal & ex = dynamic_cast<const L2VectorPrimal&>(x);
563 const std::vector<Real>& xval = *ex.getVector();
564 return fem_->compute_L2_dot(xval,*vec_);
565 }
566
567 Real norm() const {
568 Real val = 0;
569 val = std::sqrt( dot(*this) );
570 return val;
571 }
572
573 ROL::Ptr<ROL::Vector<Real> > clone() const {
574 return ROL::makePtr<L2VectorPrimal>( ROL::makePtr<std::vector<Real>>(vec_->size(),0.0),fem_);
575 }
576
577 ROL::Ptr<const std::vector<Real> > getVector() const {
578 return vec_;
579 }
580
581 ROL::Ptr<std::vector<Real> > getVector() {
582 return vec_;
583 }
584
585 ROL::Ptr<ROL::Vector<Real> > basis( const int i ) const {
586 ROL::Ptr<L2VectorPrimal> e
587 = ROL::makePtr<L2VectorPrimal>( ROL::makePtr<std::vector<Real>>(vec_->size(),0.0),fem_);
588 (*e->getVector())[i] = 1.0;
589 return e;
590 }
591
592 int dimension() const {
593 return vec_->size();
594 }
595
596 const ROL::Vector<Real>& dual() const {
597 dual_vec_ = ROL::makePtr<L2VectorDual<Real>>(
598 ROL::makePtr<std::vector<Real>>(*vec_),fem_);
599
600 fem_->apply_mass(*(ROL::constPtrCast<std::vector<Real> >(dual_vec_->getVector())),*vec_);
601 return *dual_vec_;
602 }
603
604 Real apply(const ROL::Vector<Real> &x) const {
605 const L2VectorDual<Real> &ex = dynamic_cast<const L2VectorDual<Real>&>(x);
606 const std::vector<Real>& xval = *ex.getVector();
607 return std::inner_product(vec_->begin(), vec_->end(), xval.begin(), Real(0));
608 }
609
610};
611
612template<class Real>
613class L2VectorDual : public ROL::Vector<Real> {
614private:
615 ROL::Ptr<std::vector<Real> > vec_;
616 ROL::Ptr<BurgersFEM<Real> > fem_;
617
618 mutable ROL::Ptr<L2VectorPrimal<Real> > dual_vec_;
619
620public:
621 L2VectorDual(const ROL::Ptr<std::vector<Real> > & vec,
622 const ROL::Ptr<BurgersFEM<Real> > &fem)
623 : vec_(vec), fem_(fem), dual_vec_(ROL::nullPtr) {}
624
625 void set( const ROL::Vector<Real> &x ) {
626 const L2VectorDual &ex = dynamic_cast<const L2VectorDual&>(x);
627 const std::vector<Real>& xval = *ex.getVector();
628 std::copy(xval.begin(),xval.end(),vec_->begin());
629 }
630
631 void plus( const ROL::Vector<Real> &x ) {
632 const L2VectorDual &ex = dynamic_cast<const L2VectorDual&>(x);
633 const std::vector<Real>& xval = *ex.getVector();
634 unsigned dimension = vec_->size();
635 for (unsigned i=0; i<dimension; i++) {
636 (*vec_)[i] += xval[i];
637 }
638 }
639
640 void scale( const Real alpha ) {
641 unsigned dimension = vec_->size();
642 for (unsigned i=0; i<dimension; i++) {
643 (*vec_)[i] *= alpha;
644 }
645 }
646
647 Real dot( const ROL::Vector<Real> &x ) const {
648 const L2VectorDual & ex = dynamic_cast<const L2VectorDual&>(x);
649 const std::vector<Real>& xval = *ex.getVector();
650 unsigned dimension = vec_->size();
651 std::vector<Real> Mx(dimension,0.0);
652 fem_->apply_inverse_mass(Mx,xval);
653 Real val = 0.0;
654 for (unsigned i = 0; i < dimension; i++) {
655 val += Mx[i]*(*vec_)[i];
656 }
657 return val;
658 }
659
660 Real norm() const {
661 Real val = 0;
662 val = std::sqrt( dot(*this) );
663 return val;
664 }
665
666 ROL::Ptr<ROL::Vector<Real> > clone() const {
667 return ROL::makePtr<L2VectorDual>( ROL::makePtr<std::vector<Real>>(vec_->size(),0.0),fem_);
668 }
669
670 ROL::Ptr<const std::vector<Real> > getVector() const {
671 return vec_;
672 }
673
674 ROL::Ptr<std::vector<Real> > getVector() {
675 return vec_;
676 }
677
678 ROL::Ptr<ROL::Vector<Real> > basis( const int i ) const {
679 ROL::Ptr<L2VectorDual> e
680 = ROL::makePtr<L2VectorDual>( ROL::makePtr<std::vector<Real>>(vec_->size(),0.0),fem_);
681 (*e->getVector())[i] = 1.0;
682 return e;
683 }
684
685 int dimension() const {
686 return vec_->size();
687 }
688
689 const ROL::Vector<Real>& dual() const {
690 dual_vec_ = ROL::makePtr<L2VectorPrimal<Real>>(
691 ROL::makePtr<std::vector<Real>>(*vec_),fem_);
692
693 fem_->apply_inverse_mass(*(ROL::constPtrCast<std::vector<Real> >(dual_vec_->getVector())),*vec_);
694 return *dual_vec_;
695 }
696
697 Real apply(const ROL::Vector<Real> &x) const {
698 const L2VectorPrimal<Real> &ex = dynamic_cast<const L2VectorPrimal<Real>&>(x);
699 const std::vector<Real>& xval = *ex.getVector();
700 return std::inner_product(vec_->begin(), vec_->end(), xval.begin(), Real(0));
701 }
702
703};
704
705template<class Real>
706class H1VectorPrimal : public ROL::Vector<Real> {
707private:
708 ROL::Ptr<std::vector<Real> > vec_;
709 ROL::Ptr<BurgersFEM<Real> > fem_;
710
711 mutable ROL::Ptr<H1VectorDual<Real> > dual_vec_;
712
713public:
714 H1VectorPrimal(const ROL::Ptr<std::vector<Real> > & vec,
715 const ROL::Ptr<BurgersFEM<Real> > &fem)
716 : vec_(vec), fem_(fem), dual_vec_(ROL::nullPtr) {}
717
718 void set( const ROL::Vector<Real> &x ) {
719 const H1VectorPrimal &ex = dynamic_cast<const H1VectorPrimal&>(x);
720 const std::vector<Real>& xval = *ex.getVector();
721 std::copy(xval.begin(),xval.end(),vec_->begin());
722 }
723
724 void plus( const ROL::Vector<Real> &x ) {
725 const H1VectorPrimal &ex = dynamic_cast<const H1VectorPrimal&>(x);
726 const std::vector<Real>& xval = *ex.getVector();
727 unsigned dimension = vec_->size();
728 for (unsigned i=0; i<dimension; i++) {
729 (*vec_)[i] += xval[i];
730 }
731 }
732
733 void scale( const Real alpha ) {
734 unsigned dimension = vec_->size();
735 for (unsigned i=0; i<dimension; i++) {
736 (*vec_)[i] *= alpha;
737 }
738 }
739
740 Real dot( const ROL::Vector<Real> &x ) const {
741 const H1VectorPrimal & ex = dynamic_cast<const H1VectorPrimal&>(x);
742 const std::vector<Real>& xval = *ex.getVector();
743 return fem_->compute_H1_dot(xval,*vec_);
744 }
745
746 Real norm() const {
747 Real val = 0;
748 val = std::sqrt( dot(*this) );
749 return val;
750 }
751
752 ROL::Ptr<ROL::Vector<Real> > clone() const {
753 return ROL::makePtr<H1VectorPrimal>( ROL::makePtr<std::vector<Real>>(vec_->size(),0.0),fem_);
754 }
755
756 ROL::Ptr<const std::vector<Real> > getVector() const {
757 return vec_;
758 }
759
760 ROL::Ptr<std::vector<Real> > getVector() {
761 return vec_;
762 }
763
764 ROL::Ptr<ROL::Vector<Real> > basis( const int i ) const {
765 ROL::Ptr<H1VectorPrimal> e
766 = ROL::makePtr<H1VectorPrimal>( ROL::makePtr<std::vector<Real>>(vec_->size(),0.0),fem_);
767 (*e->getVector())[i] = 1.0;
768 return e;
769 }
770
771 int dimension() const {
772 return vec_->size();
773 }
774
775 const ROL::Vector<Real>& dual() const {
776 dual_vec_ = ROL::makePtr<H1VectorDual<Real>>(
777 ROL::makePtr<std::vector<Real>>(*vec_),fem_);
778
779 fem_->apply_H1(*(ROL::constPtrCast<std::vector<Real> >(dual_vec_->getVector())),*vec_);
780 return *dual_vec_;
781 }
782
783 Real apply(const ROL::Vector<Real> &x) const {
784 const H1VectorDual<Real> &ex = dynamic_cast<const H1VectorDual<Real>&>(x);
785 const std::vector<Real>& xval = *ex.getVector();
786 return std::inner_product(vec_->begin(), vec_->end(), xval.begin(), Real(0));
787 }
788
789};
790
791template<class Real>
792class H1VectorDual : public ROL::Vector<Real> {
793private:
794 ROL::Ptr<std::vector<Real> > vec_;
795 ROL::Ptr<BurgersFEM<Real> > fem_;
796
797 mutable ROL::Ptr<H1VectorPrimal<Real> > dual_vec_;
798
799public:
800 H1VectorDual(const ROL::Ptr<std::vector<Real> > & vec,
801 const ROL::Ptr<BurgersFEM<Real> > &fem)
802 : vec_(vec), fem_(fem), dual_vec_(ROL::nullPtr) {}
803
804 void set( const ROL::Vector<Real> &x ) {
805 const H1VectorDual &ex = dynamic_cast<const H1VectorDual&>(x);
806 const std::vector<Real>& xval = *ex.getVector();
807 std::copy(xval.begin(),xval.end(),vec_->begin());
808 }
809
810 void plus( const ROL::Vector<Real> &x ) {
811 const H1VectorDual &ex = dynamic_cast<const H1VectorDual&>(x);
812 const std::vector<Real>& xval = *ex.getVector();
813 unsigned dimension = vec_->size();
814 for (unsigned i=0; i<dimension; i++) {
815 (*vec_)[i] += xval[i];
816 }
817 }
818
819 void scale( const Real alpha ) {
820 unsigned dimension = vec_->size();
821 for (unsigned i=0; i<dimension; i++) {
822 (*vec_)[i] *= alpha;
823 }
824 }
825
826 Real dot( const ROL::Vector<Real> &x ) const {
827 const H1VectorDual & ex = dynamic_cast<const H1VectorDual&>(x);
828 const std::vector<Real>& xval = *ex.getVector();
829 unsigned dimension = vec_->size();
830 std::vector<Real> Mx(dimension,0.0);
831 fem_->apply_inverse_H1(Mx,xval);
832 Real val = 0.0;
833 for (unsigned i = 0; i < dimension; i++) {
834 val += Mx[i]*(*vec_)[i];
835 }
836 return val;
837 }
838
839 Real norm() const {
840 Real val = 0;
841 val = std::sqrt( dot(*this) );
842 return val;
843 }
844
845 ROL::Ptr<ROL::Vector<Real> > clone() const {
846 return ROL::makePtr<H1VectorDual>( ROL::makePtr<std::vector<Real>>(vec_->size(),0.0),fem_);
847 }
848
849 ROL::Ptr<const std::vector<Real> > getVector() const {
850 return vec_;
851 }
852
853 ROL::Ptr<std::vector<Real> > getVector() {
854 return vec_;
855 }
856
857 ROL::Ptr<ROL::Vector<Real> > basis( const int i ) const {
858 ROL::Ptr<H1VectorDual> e
859 = ROL::makePtr<H1VectorDual>( ROL::makePtr<std::vector<Real>>(vec_->size(),0.0),fem_);
860 (*e->getVector())[i] = 1.0;
861 return e;
862 }
863
864 int dimension() const {
865 return vec_->size();
866 }
867
868 const ROL::Vector<Real>& dual() const {
869 dual_vec_ = ROL::makePtr<H1VectorPrimal<Real>>(
870 ROL::makePtr<std::vector<Real>>(*vec_),fem_);
871
872 fem_->apply_inverse_H1(*(ROL::constPtrCast<std::vector<Real> >(dual_vec_->getVector())),*vec_);
873 return *dual_vec_;
874 }
875
876 Real apply(const ROL::Vector<Real> &x) const {
877 const H1VectorPrimal<Real> &ex = dynamic_cast<const H1VectorPrimal<Real>&>(x);
878 const std::vector<Real>& xval = *ex.getVector();
879 return std::inner_product(vec_->begin(), vec_->end(), xval.begin(), Real(0));
880 }
881
882};
883
884template<class Real>
886private:
887
890
893
896
897 ROL::Ptr<BurgersFEM<Real> > fem_;
898 bool useHessian_;
899
900public:
901 Constraint_BurgersControl(ROL::Ptr<BurgersFEM<Real> > &fem, bool useHessian = true)
902 : fem_(fem), useHessian_(useHessian) {}
903
905 const ROL::Vector<Real> &z, Real &tol) {
906 ROL::Ptr<std::vector<Real> > cp =
907 dynamic_cast<PrimalConstraintVector&>(c).getVector();
908 ROL::Ptr<const std::vector<Real> > up =
909 dynamic_cast<const PrimalStateVector&>(u).getVector();
910 ROL::Ptr<const std::vector<Real> > zp =
911 dynamic_cast<const PrimalControlVector&>(z).getVector();
912 fem_->compute_residual(*cp,*up,*zp);
913 }
914
916
918 const ROL::Vector<Real> &z, Real &tol) {
919 ROL::Ptr<std::vector<Real> > jvp =
920 dynamic_cast<PrimalConstraintVector&>(jv).getVector();
921 ROL::Ptr<const std::vector<Real> > vp =
922 dynamic_cast<const PrimalStateVector&>(v).getVector();
923 ROL::Ptr<const std::vector<Real> > up =
924 dynamic_cast<const PrimalStateVector&>(u).getVector();
925 ROL::Ptr<const std::vector<Real> > zp =
926 dynamic_cast<const PrimalControlVector&>(z).getVector();
927 fem_->apply_pde_jacobian(*jvp,*vp,*up,*zp);
928 }
929
931 const ROL::Vector<Real> &z, Real &tol) {
932 ROL::Ptr<std::vector<Real> > jvp =
933 dynamic_cast<PrimalConstraintVector&>(jv).getVector();
934 ROL::Ptr<const std::vector<Real> > vp =
935 dynamic_cast<const PrimalControlVector&>(v).getVector();
936 ROL::Ptr<const std::vector<Real> > up =
937 dynamic_cast<const PrimalStateVector&>(u).getVector();
938 ROL::Ptr<const std::vector<Real> > zp =
939 dynamic_cast<const PrimalControlVector&>(z).getVector();
940 fem_->apply_control_jacobian(*jvp,*vp,*up,*zp);
941 }
942
944 const ROL::Vector<Real> &z, Real &tol) {
945 ROL::Ptr<std::vector<Real> > ijvp =
946 dynamic_cast<PrimalStateVector&>(ijv).getVector();
947 ROL::Ptr<const std::vector<Real> > vp =
948 dynamic_cast<const PrimalConstraintVector&>(v).getVector();
949 ROL::Ptr<const std::vector<Real> > up =
950 dynamic_cast<const PrimalStateVector&>(u).getVector();
951 ROL::Ptr<const std::vector<Real> > zp =
952 dynamic_cast<const PrimalControlVector&>(z).getVector();
953 fem_->apply_inverse_pde_jacobian(*ijvp,*vp,*up,*zp);
954 }
955
957 const ROL::Vector<Real> &z, Real &tol) {
958 ROL::Ptr<std::vector<Real> > jvp =
959 dynamic_cast<DualStateVector&>(ajv).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 fem_->apply_adjoint_pde_jacobian(*jvp,*vp,*up,*zp);
967 }
968
970 const ROL::Vector<Real> &z, Real &tol) {
971 ROL::Ptr<std::vector<Real> > jvp =
972 dynamic_cast<DualControlVector&>(jv).getVector();
973 ROL::Ptr<const std::vector<Real> > vp =
974 dynamic_cast<const DualConstraintVector&>(v).getVector();
975 ROL::Ptr<const std::vector<Real> > up =
976 dynamic_cast<const PrimalStateVector&>(u).getVector();
977 ROL::Ptr<const std::vector<Real> > zp =
978 dynamic_cast<const PrimalControlVector&>(z).getVector();
979 fem_->apply_adjoint_control_jacobian(*jvp,*vp,*up,*zp);
980 }
981
983 const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
984 ROL::Ptr<std::vector<Real> > iajvp =
985 dynamic_cast<DualConstraintVector&>(iajv).getVector();
986 ROL::Ptr<const std::vector<Real> > vp =
987 dynamic_cast<const DualStateVector&>(v).getVector();
988 ROL::Ptr<const std::vector<Real> > up =
989 dynamic_cast<const PrimalStateVector&>(u).getVector();
990 ROL::Ptr<const std::vector<Real> > zp =
991 dynamic_cast<const PrimalControlVector&>(z).getVector();
992 fem_->apply_inverse_adjoint_pde_jacobian(*iajvp,*vp,*up,*zp);
993 }
994
996 const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
997 if ( useHessian_ ) {
998 ROL::Ptr<std::vector<Real> > ahwvp =
999 dynamic_cast<DualStateVector&>(ahwv).getVector();
1000 ROL::Ptr<const std::vector<Real> > wp =
1001 dynamic_cast<const DualConstraintVector&>(w).getVector();
1002 ROL::Ptr<const std::vector<Real> > vp =
1003 dynamic_cast<const PrimalStateVector&>(v).getVector();
1004 ROL::Ptr<const std::vector<Real> > up =
1005 dynamic_cast<const PrimalStateVector&>(u).getVector();
1006 ROL::Ptr<const std::vector<Real> > zp =
1007 dynamic_cast<const PrimalControlVector&>(z).getVector();
1008 fem_->apply_adjoint_pde_hessian(*ahwvp,*wp,*vp,*up,*zp);
1009 }
1010 else {
1011 ahwv.zero();
1012 }
1013 }
1014
1016 const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
1017 if ( useHessian_ ) {
1018 ROL::Ptr<std::vector<Real> > ahwvp =
1019 dynamic_cast<DualControlVector&>(ahwv).getVector();
1020 ROL::Ptr<const std::vector<Real> > wp =
1021 dynamic_cast<const DualConstraintVector&>(w).getVector();
1022 ROL::Ptr<const std::vector<Real> > vp =
1023 dynamic_cast<const PrimalStateVector&>(v).getVector();
1024 ROL::Ptr<const std::vector<Real> > up =
1025 dynamic_cast<const PrimalStateVector&>(u).getVector();
1026 ROL::Ptr<const std::vector<Real> > zp =
1027 dynamic_cast<const PrimalControlVector&>(z).getVector();
1028 fem_->apply_adjoint_control_pde_hessian(*ahwvp,*wp,*vp,*up,*zp);
1029 }
1030 else {
1031 ahwv.zero();
1032 }
1033 }
1035 const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
1036 if ( useHessian_ ) {
1037 ROL::Ptr<std::vector<Real> > ahwvp =
1038 dynamic_cast<DualStateVector&>(ahwv).getVector();
1039 ROL::Ptr<const std::vector<Real> > wp =
1040 dynamic_cast<const DualConstraintVector&>(w).getVector();
1041 ROL::Ptr<const std::vector<Real> > vp =
1042 dynamic_cast<const PrimalControlVector&>(v).getVector();
1043 ROL::Ptr<const std::vector<Real> > up =
1044 dynamic_cast<const PrimalStateVector&>(u).getVector();
1045 ROL::Ptr<const std::vector<Real> > zp =
1046 dynamic_cast<const PrimalControlVector&>(z).getVector();
1047 fem_->apply_adjoint_pde_control_hessian(*ahwvp,*wp,*vp,*up,*zp);
1048 }
1049 else {
1050 ahwv.zero();
1051 }
1052 }
1054 const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
1055 if ( useHessian_ ) {
1056 ROL::Ptr<std::vector<Real> > ahwvp =
1057 dynamic_cast<DualControlVector&>(ahwv).getVector();
1058 ROL::Ptr<const std::vector<Real> > wp =
1059 dynamic_cast<const DualConstraintVector&>(w).getVector();
1060 ROL::Ptr<const std::vector<Real> > vp =
1061 dynamic_cast<const PrimalControlVector&>(v).getVector();
1062 ROL::Ptr<const std::vector<Real> > up =
1063 dynamic_cast<const PrimalStateVector&>(u).getVector();
1064 ROL::Ptr<const std::vector<Real> > zp =
1065 dynamic_cast<const PrimalControlVector&>(z).getVector();
1066 fem_->apply_adjoint_control_hessian(*ahwvp,*wp,*vp,*up,*zp);
1067 }
1068 else {
1069 ahwv.zero();
1070 }
1071 }
1072};
1073
1074template<class Real>
1076private:
1077
1080
1083
1084 Real alpha_; // Penalty Parameter
1085 ROL::Ptr<BurgersFEM<Real> > fem_;
1086 ROL::Ptr<ROL::Vector<Real> > ud_;
1087 ROL::Ptr<ROL::Vector<Real> > diff_;
1088
1089public:
1091 const ROL::Ptr<ROL::Vector<Real> > &ud,
1092 Real alpha = 1.e-4) : alpha_(alpha), fem_(fem), ud_(ud) {
1093 diff_ = ud_->clone();
1094 }
1095
1096 Real value( const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
1097 ROL::Ptr<const std::vector<Real> > up =
1098 dynamic_cast<const PrimalStateVector&>(u).getVector();
1099 ROL::Ptr<const std::vector<Real> > zp =
1100 dynamic_cast<const PrimalControlVector&>(z).getVector();
1101 ROL::Ptr<const std::vector<Real> > udp =
1102 dynamic_cast<const L2VectorPrimal<Real>&>(*ud_).getVector();
1103
1104 std::vector<Real> diff(udp->size(),0.0);
1105 for (unsigned i = 0; i < udp->size(); i++) {
1106 diff[i] = (*up)[i] - (*udp)[i];
1107 }
1108 return 0.5*(fem_->compute_L2_dot(diff,diff) + alpha_*fem_->compute_L2_dot(*zp,*zp));
1109 }
1110
1111 void gradient_1( ROL::Vector<Real> &g, const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
1112 ROL::Ptr<std::vector<Real> > gp =
1113 dynamic_cast<DualStateVector&>(g).getVector();
1114 ROL::Ptr<const std::vector<Real> > up =
1115 dynamic_cast<const PrimalStateVector&>(u).getVector();
1116 ROL::Ptr<const std::vector<Real> > udp =
1117 dynamic_cast<const L2VectorPrimal<Real>&>(*ud_).getVector();
1118
1119 std::vector<Real> diff(udp->size(),0.0);
1120 for (unsigned i = 0; i < udp->size(); i++) {
1121 diff[i] = (*up)[i] - (*udp)[i];
1122 }
1123 fem_->apply_mass(*gp,diff);
1124 }
1125
1126 void gradient_2( ROL::Vector<Real> &g, const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
1127 ROL::Ptr<std::vector<Real> > gp =
1128 dynamic_cast<DualControlVector&>(g).getVector();
1129 ROL::Ptr<const std::vector<Real> > zp =
1130 dynamic_cast<const PrimalControlVector&>(z).getVector();
1131
1132 fem_->apply_mass(*gp,*zp);
1133 for (unsigned i = 0; i < zp->size(); i++) {
1134 (*gp)[i] *= alpha_;
1135 }
1136 }
1137
1139 const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
1140 ROL::Ptr<std::vector<Real> > hvp =
1141 dynamic_cast<DualStateVector&>(hv).getVector();
1142 ROL::Ptr<const std::vector<Real> > vp =
1143 dynamic_cast<const PrimalStateVector&>(v).getVector();
1144
1145 fem_->apply_mass(*hvp,*vp);
1146 }
1147
1149 const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
1150 hv.zero();
1151 }
1152
1154 const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
1155 hv.zero();
1156 }
1157
1159 const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
1160 ROL::Ptr<std::vector<Real> > hvp =
1161 dynamic_cast<DualControlVector&>(hv).getVector();
1162 ROL::Ptr<const std::vector<Real> > vp =
1163 dynamic_cast<const PrimalControlVector&>(v).getVector();
1164
1165 fem_->apply_mass(*hvp,*vp);
1166 for (unsigned i = 0; i < vp->size(); i++) {
1167 (*hvp)[i] *= alpha_;
1168 }
1169 }
1170};
1171
1172template<class Real>
1174private:
1175 int dim_;
1176 std::vector<Real> x_lo_;
1177 std::vector<Real> x_up_;
1180 ROL::Ptr<BurgersFEM<Real> > fem_;
1181 ROL::Ptr<ROL::Vector<Real> > l_;
1182 ROL::Ptr<ROL::Vector<Real> > u_;
1183
1184 void cast_vector(ROL::Ptr<std::vector<Real> > &xvec,
1185 ROL::Vector<Real> &x) const {
1186 try {
1187 xvec = dynamic_cast<L2VectorPrimal<Real>&>(x).getVector();
1188 }
1189 catch (std::exception &e) {
1190 xvec = dynamic_cast<L2VectorDual<Real>&>(x).getVector();
1191 }
1192 }
1193
1194 void cast_const_vector(ROL::Ptr<const std::vector<Real> > &xvec,
1195 const ROL::Vector<Real> &x) const {
1196 try {
1197 xvec = dynamic_cast<const L2VectorPrimal<Real>&>(x).getVector();
1198 }
1199 catch (std::exception &e) {
1200 xvec = dynamic_cast<const L2VectorDual<Real>&>(x).getVector();
1201 }
1202 }
1203
1204 void axpy(std::vector<Real> &out, const Real a,
1205 const std::vector<Real> &x, const std::vector<Real> &y) const{
1206 out.resize(dim_,0.0);
1207 for (unsigned i = 0; i < dim_; i++) {
1208 out[i] = a*x[i] + y[i];
1209 }
1210 }
1211
1212 void projection(std::vector<Real> &x) {
1213 for ( int i = 0; i < dim_; i++ ) {
1214 x[i] = std::max(x_lo_[i],std::min(x_up_[i],x[i]));
1215 }
1216 }
1217
1218public:
1219 L2BoundConstraint(std::vector<Real> &l, std::vector<Real> &u,
1220 const ROL::Ptr<BurgersFEM<Real> > &fem, Real scale = 1.0)
1221 : x_lo_(l), x_up_(u), scale_(scale), fem_(fem) {
1222 dim_ = x_lo_.size();
1223 for ( int i = 0; i < dim_; i++ ) {
1224 if ( i == 0 ) {
1225 min_diff_ = x_up_[i] - x_lo_[i];
1226 }
1227 else {
1228 min_diff_ = ( (min_diff_ < (x_up_[i] - x_lo_[i])) ? min_diff_ : (x_up_[i] - x_lo_[i]) );
1229 }
1230 }
1231 min_diff_ *= 0.5;
1232 l_ = ROL::makePtr<L2VectorPrimal<Real>>(
1233 ROL::makePtr<std::vector<Real>>(l), fem);
1234 u_ = ROL::makePtr<L2VectorPrimal<Real>>(
1235 ROL::makePtr<std::vector<Real>>(u), fem);
1236 }
1237
1239 ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1240 bool val = true;
1241 int cnt = 1;
1242 for ( int i = 0; i < dim_; i++ ) {
1243 if ( (*ex)[i] >= x_lo_[i] && (*ex)[i] <= x_up_[i] ) { cnt *= 1; }
1244 else { cnt *= 0; }
1245 }
1246 if ( cnt == 0 ) { val = false; }
1247 return val;
1248 }
1249
1251 ROL::Ptr<std::vector<Real> > ex; cast_vector(ex,x);
1252 projection(*ex);
1253 }
1254
1255 void pruneLowerActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &x, Real eps = Real(0)) {
1256 ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1257 ROL::Ptr<std::vector<Real> > ev; cast_vector(ev,v);
1258 Real epsn = std::min(scale_*eps,min_diff_);
1259 for ( int i = 0; i < dim_; i++ ) {
1260 if ( ((*ex)[i] <= x_lo_[i]+epsn) ) {
1261 (*ev)[i] = 0.0;
1262 }
1263 }
1264 }
1265
1266 void pruneUpperActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &x, Real eps = Real(0)) {
1267 ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1268 ROL::Ptr<std::vector<Real> > ev; cast_vector(ev,v);
1269 Real epsn = std::min(scale_*eps,min_diff_);
1270 for ( int i = 0; i < dim_; i++ ) {
1271 if ( ((*ex)[i] >= x_up_[i]-epsn) ) {
1272 (*ev)[i] = 0.0;
1273 }
1274 }
1275 }
1276
1277 void pruneActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &x, Real eps = Real(0)) {
1278 ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1279 ROL::Ptr<std::vector<Real> > ev; cast_vector(ev,v);
1280 Real epsn = std::min(scale_*eps,min_diff_);
1281 for ( int i = 0; i < dim_; i++ ) {
1282 if ( ((*ex)[i] <= x_lo_[i]+epsn) ||
1283 ((*ex)[i] >= x_up_[i]-epsn) ) {
1284 (*ev)[i] = 0.0;
1285 }
1286 }
1287 }
1288
1289 void pruneLowerActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &g, const ROL::Vector<Real> &x, Real xeps = Real(0), Real geps = Real(0)) {
1290 ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1291 ROL::Ptr<const std::vector<Real> > eg; cast_const_vector(eg,g);
1292 ROL::Ptr<std::vector<Real> > ev; cast_vector(ev,v);
1293 Real epsn = std::min(scale_*xeps,min_diff_);
1294 for ( int i = 0; i < dim_; i++ ) {
1295 if ( ((*ex)[i] <= x_lo_[i]+epsn && (*eg)[i] > geps) ) {
1296 (*ev)[i] = 0.0;
1297 }
1298 }
1299 }
1300
1301 void pruneUpperActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &g, const ROL::Vector<Real> &x, Real xeps = Real(0), Real geps = Real(0)) {
1302 ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1303 ROL::Ptr<const std::vector<Real> > eg; cast_const_vector(eg,g);
1304 ROL::Ptr<std::vector<Real> > ev; cast_vector(ev,v);
1305 Real epsn = std::min(scale_*xeps,min_diff_);
1306 for ( int i = 0; i < dim_; i++ ) {
1307 if ( ((*ex)[i] >= x_up_[i]-epsn && (*eg)[i] < geps) ) {
1308 (*ev)[i] = 0.0;
1309 }
1310 }
1311 }
1312
1313 void pruneActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &g, const ROL::Vector<Real> &x, Real xeps = Real(0), Real geps = Real(0)) {
1314 ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1315 ROL::Ptr<const std::vector<Real> > eg; cast_const_vector(eg,g);
1316 ROL::Ptr<std::vector<Real> > ev; cast_vector(ev,v);
1317 Real epsn = std::min(scale_*xeps,min_diff_);
1318 for ( int i = 0; i < dim_; i++ ) {
1319 if ( ((*ex)[i] <= x_lo_[i]+epsn && (*eg)[i] > geps) ||
1320 ((*ex)[i] >= x_up_[i]-epsn && (*eg)[i] < -geps) ) {
1321 (*ev)[i] = 0.0;
1322 }
1323 }
1324 }
1325
1326 const ROL::Ptr<const ROL::Vector<Real> > getLowerBound(void) const {
1327 return l_;
1328 }
1329
1330 const ROL::Ptr<const ROL::Vector<Real> > getUpperBound(void) const {
1331 return u_;
1332 }
1333};
1334
1335template<class Real>
1337private:
1338 int dim_;
1339 std::vector<Real> x_lo_;
1340 std::vector<Real> x_up_;
1343 ROL::Ptr<BurgersFEM<Real> > fem_;
1344 ROL::Ptr<ROL::Vector<Real> > l_;
1345 ROL::Ptr<ROL::Vector<Real> > u_;
1346
1347 void cast_vector(ROL::Ptr<std::vector<Real> > &xvec,
1348 ROL::Vector<Real> &x) const {
1349 try {
1350 xvec = dynamic_cast<H1VectorPrimal<Real>&>(x).getVector();
1351 }
1352 catch (std::exception &e) {
1353 xvec = dynamic_cast<H1VectorDual<Real>&>(x).getVector();
1354 }
1355 }
1356
1357 void cast_const_vector(ROL::Ptr<const std::vector<Real> > &xvec,
1358 const ROL::Vector<Real> &x) const {
1359 try {
1360 xvec = dynamic_cast<const H1VectorPrimal<Real>&>(x).getVector();
1361 }
1362 catch (std::exception &e) {
1363 xvec = dynamic_cast<const H1VectorDual<Real>&>(x).getVector();
1364 }
1365 }
1366
1367 void axpy(std::vector<Real> &out, const Real a,
1368 const std::vector<Real> &x, const std::vector<Real> &y) const{
1369 out.resize(dim_,0.0);
1370 for (unsigned i = 0; i < dim_; i++) {
1371 out[i] = a*x[i] + y[i];
1372 }
1373 }
1374
1375 void projection(std::vector<Real> &x) {
1376 for ( int i = 0; i < dim_; i++ ) {
1377 x[i] = std::max(x_lo_[i],std::min(x_up_[i],x[i]));
1378 }
1379 }
1380
1381public:
1382 H1BoundConstraint(std::vector<Real> &l, std::vector<Real> &u,
1383 const ROL::Ptr<BurgersFEM<Real> > &fem, Real scale = 1.0)
1384 : x_lo_(l), x_up_(u), scale_(scale), fem_(fem) {
1385 dim_ = x_lo_.size();
1386 for ( int i = 0; i < dim_; i++ ) {
1387 if ( i == 0 ) {
1388 min_diff_ = x_up_[i] - x_lo_[i];
1389 }
1390 else {
1391 min_diff_ = ( (min_diff_ < (x_up_[i] - x_lo_[i])) ? min_diff_ : (x_up_[i] - x_lo_[i]) );
1392 }
1393 }
1394 min_diff_ *= 0.5;
1395 l_ = ROL::makePtr<H1VectorPrimal<Real>>(
1396 ROL::makePtr<std::vector<Real>>(l), fem);
1397 u_ = ROL::makePtr<H1VectorPrimal<Real>>(
1398 ROL::makePtr<std::vector<Real>>(u), fem);
1399 }
1400
1402 ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1403 bool val = true;
1404 int cnt = 1;
1405 for ( int i = 0; i < dim_; i++ ) {
1406 if ( (*ex)[i] >= x_lo_[i] && (*ex)[i] <= x_up_[i] ) { cnt *= 1; }
1407 else { cnt *= 0; }
1408 }
1409 if ( cnt == 0 ) { val = false; }
1410 return val;
1411 }
1412
1414 ROL::Ptr<std::vector<Real> > ex; cast_vector(ex,x);
1415 projection(*ex);
1416 }
1417
1418 void pruneLowerActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &x, Real eps = Real(0)) {
1419 ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1420 ROL::Ptr<std::vector<Real> > ev; cast_vector(ev,v);
1421 Real epsn = std::min(scale_*eps,min_diff_);
1422 for ( int i = 0; i < dim_; i++ ) {
1423 if ( ((*ex)[i] <= x_lo_[i]+epsn) ) {
1424 (*ev)[i] = 0.0;
1425 }
1426 }
1427 }
1428
1429 void pruneUpperActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &x, Real eps = Real(0)) {
1430 ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1431 ROL::Ptr<std::vector<Real> > ev; cast_vector(ev,v);
1432 Real epsn = std::min(scale_*eps,min_diff_);
1433 for ( int i = 0; i < dim_; i++ ) {
1434 if ( ((*ex)[i] >= x_up_[i]-epsn) ) {
1435 (*ev)[i] = 0.0;
1436 }
1437 }
1438 }
1439
1440 void pruneActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &x, Real eps = Real(0)) {
1441 ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1442 ROL::Ptr<std::vector<Real> > ev; cast_vector(ev,v);
1443 Real epsn = std::min(scale_*eps,min_diff_);
1444 for ( int i = 0; i < dim_; i++ ) {
1445 if ( ((*ex)[i] <= x_lo_[i]+epsn) ||
1446 ((*ex)[i] >= x_up_[i]-epsn) ) {
1447 (*ev)[i] = 0.0;
1448 }
1449 }
1450 }
1451
1452 void pruneLowerActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &g, const ROL::Vector<Real> &x, Real xeps = Real(0), Real geps = Real(0)) {
1453 ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1454 ROL::Ptr<const std::vector<Real> > eg; cast_const_vector(eg,g);
1455 ROL::Ptr<std::vector<Real> > ev; cast_vector(ev,v);
1456 Real epsn = std::min(scale_*xeps,min_diff_);
1457 for ( int i = 0; i < dim_; i++ ) {
1458 if ( ((*ex)[i] <= x_lo_[i]+epsn && (*eg)[i] > geps) ) {
1459 (*ev)[i] = 0.0;
1460 }
1461 }
1462 }
1463
1464 void pruneUpperActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &g, const ROL::Vector<Real> &x, Real xeps = Real(0), Real geps = Real(0)) {
1465 ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1466 ROL::Ptr<const std::vector<Real> > eg; cast_const_vector(eg,g);
1467 ROL::Ptr<std::vector<Real> > ev; cast_vector(ev,v);
1468 Real epsn = std::min(scale_*xeps,min_diff_);
1469 for ( int i = 0; i < dim_; i++ ) {
1470 if ( ((*ex)[i] >= x_up_[i]-epsn && (*eg)[i] < -geps) ) {
1471 (*ev)[i] = 0.0;
1472 }
1473 }
1474 }
1475
1476 void pruneActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &g, const ROL::Vector<Real> &x, Real xeps = Real(0), Real geps = Real(0)) {
1477 ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1478 ROL::Ptr<const std::vector<Real> > eg; cast_const_vector(eg,g);
1479 ROL::Ptr<std::vector<Real> > ev; cast_vector(ev,v);
1480 Real epsn = std::min(scale_*xeps,min_diff_);
1481 for ( int i = 0; i < dim_; i++ ) {
1482 if ( ((*ex)[i] <= x_lo_[i]+epsn && (*eg)[i] > geps) ||
1483 ((*ex)[i] >= x_up_[i]-epsn && (*eg)[i] < -geps) ) {
1484 (*ev)[i] = 0.0;
1485 }
1486 }
1487 }
1488
1489 const ROL::Ptr<const ROL::Vector<Real> > getLowerBound(void) const {
1490 return l_;
1491 }
1492
1493 const ROL::Ptr<const ROL::Vector<Real> > getUpperBound(void) const {
1494 return u_;
1495 }
1496};
Contains definitions of custom data types in ROL.
BurgersFEM(int nx=128, Real nu=1.e-2, Real nl=1.0, Real u0=1.0, Real u1=0.0, Real f=0.0, Real cH1=1.0, Real cL2=1.0)
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
void apply_inverse_H1(std::vector< Real > &Mu, const std::vector< Real > &u) const
Real compute_H1_norm(const std::vector< Real > &r) const
Real compute_L2_dot(const std::vector< Real > &x, const std::vector< Real > &y) const
void scale(std::vector< Real > &u, const Real alpha=0.0) const
void apply_inverse_mass(std::vector< Real > &Mu, const std::vector< Real > &u) const
void compute_residual(std::vector< Real > &r, const std::vector< Real > &u, const std::vector< Real > &z) const
void apply_control_jacobian(std::vector< Real > &jv, const std::vector< Real > &v, const std::vector< Real > &u, const std::vector< Real > &z) const
void axpy(std::vector< Real > &out, const Real a, const std::vector< Real > &x, const std::vector< Real > &y) const
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
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)
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)
Real mesh_spacing(void) const
void test_inverse_mass(std::ostream &outStream=std::cout)
void test_inverse_H1(std::ostream &outStream=std::cout)
void apply_mass(std::vector< Real > &Mu, const std::vector< Real > &u) const
void apply_H1(std::vector< Real > &Mu, const std::vector< Real > &u) const
int num_dof(void) const
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
Real compute_H1_dot(const std::vector< Real > &x, const std::vector< Real > &y) const
void apply_pde_jacobian(std::vector< Real > &jv, const std::vector< Real > &v, const std::vector< Real > &u, const std::vector< Real > &z) const
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
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)
void compute_pde_jacobian(std::vector< Real > &dl, std::vector< Real > &d, std::vector< Real > &du, const std::vector< Real > &u) const
Real compute_L2_norm(const std::vector< Real > &r) const
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
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
Real cH1_
Definition test_04.hpp:78
void update(std::vector< Real > &u, const std::vector< Real > &s, const Real alpha=1.0) const
Real cL2_
Definition test_04.hpp:79
L2VectorPrimal< Real > PrimalControlVector
H1VectorDual< Real > PrimalConstraintVector
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 ...
void applyAdjointHessian_22(ROL::Vector< Real > &ahwv, const ROL::Vector< Real > &w, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Apply the optimization-space derivative of the adjoint of the constraint optimization-space Jacobian ...
void 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...
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...
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 .
ROL::Ptr< BurgersFEM< Real > > fem_
Definition test_04.hpp:716
H1VectorPrimal< Real > PrimalStateVector
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...
H1VectorDual< Real > DualStateVector
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 .
H1VectorPrimal< Real > DualConstraintVector
L2VectorDual< Real > DualControlVector
void value(ROL::Vector< Real > &c, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Evaluate the constraint operator at .
void applyAdjointHessian_12(ROL::Vector< Real > &ahwv, const ROL::Vector< Real > &w, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Apply the optimization-space derivative of the adjoint of the constraint simulation-space Jacobian at...
void 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 .
Constraint_BurgersControl(ROL::Ptr< BurgersFEM< Real > > &fem, bool useHessian=true)
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 .
H1BoundConstraint(std::vector< Real > &l, std::vector< Real > &u, const ROL::Ptr< BurgersFEM< Real > > &fem, Real scale=1.0)
void pruneLowerActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &x, Real eps=Real(0))
Set variables to zero if they correspond to the lower -active set.
ROL::Ptr< ROL::Vector< Real > > l_
ROL::Ptr< BurgersFEM< Real > > fem_
std::vector< Real > x_lo_
void pruneUpperActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &x, Real eps=Real(0))
Set variables to zero if they correspond to the upper -active set.
void projection(std::vector< Real > &x)
std::vector< Real > x_up_
bool isFeasible(const ROL::Vector< Real > &x)
Check if the vector, v, is feasible.
void pruneUpperActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &g, const ROL::Vector< Real > &x, Real xeps=Real(0), Real geps=Real(0))
Set variables to zero if they correspond to the upper -binding set.
void axpy(std::vector< Real > &out, const Real a, const std::vector< Real > &x, const std::vector< Real > &y) const
void pruneActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &x, Real eps=Real(0))
void project(ROL::Vector< Real > &x)
Project optimization variables onto the bounds.
void pruneActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &g, const ROL::Vector< Real > &x, Real xeps=Real(0), Real geps=Real(0))
ROL::Ptr< ROL::Vector< Real > > u_
void pruneLowerActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &g, const ROL::Vector< Real > &x, Real xeps=Real(0), Real geps=Real(0))
Set variables to zero if they correspond to the -binding set.
void cast_const_vector(ROL::Ptr< const std::vector< Real > > &xvec, const ROL::Vector< Real > &x) const
const ROL::Ptr< const ROL::Vector< Real > > getLowerBound(void) const
Return the ref count pointer to the lower bound vector.
void cast_vector(ROL::Ptr< std::vector< Real > > &xvec, ROL::Vector< Real > &x) const
const ROL::Ptr< const ROL::Vector< Real > > getUpperBound(void) const
Return the ref count pointer to the upper bound vector.
Real norm() const
Returns where .
ROL::Ptr< H1VectorPrimal< Real > > dual_vec_
Definition test_04.hpp:661
int dimension() const
Return dimension of the vector space.
Real apply(const ROL::Vector< Real > &x) const
Apply to a dual vector. This is equivalent to the call .
ROL::Ptr< std::vector< Real > > vec_
Real dot(const ROL::Vector< Real > &x) const override
Compute where .
Definition test_04.hpp:670
Real dot(const ROL::Vector< Real > &x) const
Compute where .
void plus(const ROL::Vector< Real > &x)
Compute , where .
void set(const ROL::Vector< Real > &x)
Set where .
ROL::Ptr< ROL::Vector< Real > > basis(const int i) const
Return i-th basis vector.
ROL::Ptr< ROL::Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
ROL::Ptr< const std::vector< Real > > getVector() const
ROL::Ptr< BurgersFEM< Real > > fem_
Definition test_04.hpp:660
void scale(const Real alpha)
Compute where .
ROL::Ptr< std::vector< Real > > getVector()
const ROL::Vector< Real > & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis,...
H1VectorDual(const ROL::Ptr< std::vector< Real > > &vec, const ROL::Ptr< BurgersFEM< Real > > &fem)
ROL::Ptr< ROL::Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
const ROL::Vector< Real > & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis,...
void plus(const ROL::Vector< Real > &x)
Compute , where .
ROL::Ptr< const std::vector< Real > > getVector() const
void set(const ROL::Vector< Real > &x)
Set where .
ROL::Ptr< std::vector< Real > > vec_
H1VectorPrimal(const ROL::Ptr< std::vector< Real > > &vec, const ROL::Ptr< BurgersFEM< Real > > &fem)
Real norm() const
Returns where .
ROL::Ptr< ROL::Vector< Real > > basis(const int i) const
Return i-th basis vector.
ROL::Ptr< std::vector< Real > > getVector()
int dimension() const
Return dimension of the vector space.
Real dot(const ROL::Vector< Real > &x) const
Compute where .
ROL::Ptr< H1VectorDual< Real > > dual_vec_
Definition test_04.hpp:622
ROL::Ptr< BurgersFEM< Real > > fem_
Definition test_04.hpp:621
void scale(const Real alpha)
Compute where .
Real dot(const ROL::Vector< Real > &x) const override
Compute where .
Definition test_04.hpp:631
Real apply(const ROL::Vector< Real > &x) const
Apply to a dual vector. This is equivalent to the call .
void axpy(std::vector< Real > &out, const Real a, const std::vector< Real > &x, const std::vector< Real > &y) const
bool isFeasible(const ROL::Vector< Real > &x)
Check if the vector, v, is feasible.
void pruneActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &g, const ROL::Vector< Real > &x, Real xeps=Real(0), Real geps=Real(0))
ROL::Ptr< ROL::Vector< Real > > l_
L2BoundConstraint(std::vector< Real > &l, std::vector< Real > &u, const ROL::Ptr< BurgersFEM< Real > > &fem, Real scale=1.0)
void projection(std::vector< Real > &x)
const ROL::Ptr< const ROL::Vector< Real > > getLowerBound(void) const
Return the ref count pointer to the lower bound vector.
void pruneUpperActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &g, const ROL::Vector< Real > &x, Real xeps=Real(0), Real geps=Real(0))
Set variables to zero if they correspond to the upper -binding set.
void pruneActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &x, Real eps=Real(0))
std::vector< Real > x_lo_
ROL::Ptr< BurgersFEM< Real > > fem_
void pruneUpperActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &x, Real eps=Real(0))
Set variables to zero if they correspond to the upper -active set.
void pruneLowerActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &g, const ROL::Vector< Real > &x, Real xeps=Real(0), Real geps=Real(0))
Set variables to zero if they correspond to the -binding set.
std::vector< Real > x_up_
void pruneLowerActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &x, Real eps=Real(0))
Set variables to zero if they correspond to the lower -active set.
const ROL::Ptr< const ROL::Vector< Real > > getUpperBound(void) const
Return the ref count pointer to the upper bound vector.
void project(ROL::Vector< Real > &x)
Project optimization variables onto the bounds.
ROL::Ptr< ROL::Vector< Real > > u_
void cast_vector(ROL::Ptr< std::vector< Real > > &xvec, ROL::Vector< Real > &x) const
void cast_const_vector(ROL::Ptr< const std::vector< Real > > &xvec, const ROL::Vector< Real > &x) const
Real apply(const ROL::Vector< Real > &x) const
Apply to a dual vector. This is equivalent to the call .
void set(const ROL::Vector< Real > &x)
Set where .
ROL::Ptr< std::vector< Real > > getVector()
Real dot(const ROL::Vector< Real > &x) const
Compute where .
Real dot(const ROL::Vector< Real > &x) const override
Compute where .
Definition test_04.hpp:585
ROL::Ptr< ROL::Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
void scale(const Real alpha)
Compute where .
ROL::Ptr< ROL::Vector< Real > > basis(const int i) const
Return i-th basis vector.
L2VectorDual(const ROL::Ptr< std::vector< Real > > &vec, const ROL::Ptr< BurgersFEM< Real > > &fem)
ROL::Ptr< const std::vector< Real > > getVector() const
void plus(const ROL::Vector< Real > &x)
Compute , where .
int dimension() const
Return dimension of the vector space.
ROL::Ptr< std::vector< Real > > vec_
Real norm() const
Returns where .
const ROL::Vector< Real > & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis,...
ROL::Ptr< BurgersFEM< Real > > fem_
Definition test_04.hpp:575
ROL::Ptr< L2VectorPrimal< Real > > dual_vec_
Definition test_04.hpp:576
int dimension() const
Return dimension of the vector space.
ROL::Ptr< std::vector< Real > > getVector()
void set(const ROL::Vector< Real > &x)
Set where .
Real dot(const ROL::Vector< Real > &x) const
Compute where .
ROL::Ptr< const std::vector< Real > > getVector() const
const ROL::Vector< Real > & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis,...
L2VectorPrimal(const ROL::Ptr< std::vector< Real > > &vec, const ROL::Ptr< BurgersFEM< Real > > &fem)
Real norm() const
Returns where .
ROL::Ptr< ROL::Vector< Real > > basis(const int i) const
Return i-th basis vector.
ROL::Ptr< ROL::Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
void plus(const ROL::Vector< Real > &x)
Compute , where .
ROL::Ptr< L2VectorDual< Real > > dual_vec_
Definition test_04.hpp:537
ROL::Ptr< std::vector< Real > > vec_
ROL::Ptr< BurgersFEM< Real > > fem_
Definition test_04.hpp:536
Real dot(const ROL::Vector< Real > &x) const override
Compute where .
Definition test_04.hpp:546
Real apply(const ROL::Vector< Real > &x) const
Apply to a dual vector. This is equivalent to the call .
void scale(const Real alpha)
Compute where .
Objective_BurgersControl(const ROL::Ptr< BurgersFEM< Real > > &fem, const ROL::Ptr< ROL::Vector< Real > > &ud, Real alpha=1.e-4)
H1VectorPrimal< Real > PrimalStateVector
void hessVec_21(ROL::Vector< Real > &hv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
void gradient_1(ROL::Vector< Real > &g, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Compute gradient with respect to first component.
ROL::Ptr< ROL::Vector< Real > > ud_
void hessVec_11(ROL::Vector< Real > &hv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Apply Hessian approximation to vector.
void hessVec_22(ROL::Vector< Real > &hv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
ROL::Ptr< ROL::Vector< Real > > diff_
void gradient_2(ROL::Vector< Real > &g, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Compute gradient with respect to second component.
L2VectorPrimal< Real > PrimalControlVector
ROL::Ptr< BurgersFEM< Real > > fem_
H1VectorDual< Real > DualStateVector
void hessVec_12(ROL::Vector< Real > &hv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Real value(const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Compute value.
L2VectorDual< Real > DualControlVector
Provides the interface to apply upper and lower bound constraints.
Defines the constraint operator interface for simulation-based optimization.
Provides the interface to evaluate simulation-based objective functions.
Defines the linear algebra or vector space interface.
virtual void zero()
Set to zero vector.
constexpr auto dim