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