ROL
example_07.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#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 Real get_viscosity(void) const {
144 return nu_;
145 }
146
147 int num_dof(void) const {
148 return nx_;
149 }
150
151 Real mesh_spacing(void) const {
152 return dx_;
153 }
154
155 /***************************************************************************/
156 /*********************** L2 INNER PRODUCT FUNCTIONS ************************/
157 /***************************************************************************/
158 // Compute L2 inner product
159 Real compute_L2_dot(const std::vector<Real> &x, const std::vector<Real> &y) const {
160 Real ip = 0.0;
161 Real c = (((int)x.size()==nx_) ? 4.0 : 2.0);
162 for (unsigned i=0; i<x.size(); i++) {
163 if ( i == 0 ) {
164 ip += dx_/6.0*(c*x[i] + x[i+1])*y[i];
165 }
166 else if ( i == x.size()-1 ) {
167 ip += dx_/6.0*(x[i-1] + c*x[i])*y[i];
168 }
169 else {
170 ip += dx_/6.0*(x[i-1] + 4.0*x[i] + x[i+1])*y[i];
171 }
172 }
173 return ip;
174 }
175
176 // compute L2 norm
177 Real compute_L2_norm(const std::vector<Real> &r) const {
178 return std::sqrt(compute_L2_dot(r,r));
179 }
180
181 // Apply L2 Reisz operator
182 void apply_mass(std::vector<Real> &Mu, const std::vector<Real> &u ) const {
183 Mu.resize(u.size(),0.0);
184 Real c = (((int)u.size()==nx_) ? 4.0 : 2.0);
185 for (unsigned i=0; i<u.size(); i++) {
186 if ( i == 0 ) {
187 Mu[i] = dx_/6.0*(c*u[i] + u[i+1]);
188 }
189 else if ( i == u.size()-1 ) {
190 Mu[i] = dx_/6.0*(u[i-1] + c*u[i]);
191 }
192 else {
193 Mu[i] = dx_/6.0*(u[i-1] + 4.0*u[i] + u[i+1]);
194 }
195 }
196 }
197
198 // Apply L2 inverse Reisz operator
199 void apply_inverse_mass(std::vector<Real> &Mu, const std::vector<Real> &u) const {
200 unsigned nx = u.size();
201 // Build mass matrix
202 std::vector<Real> dl(nx-1,dx_/6.0);
203 std::vector<Real> d(nx,2.0*dx_/3.0);
204 std::vector<Real> du(nx-1,dx_/6.0);
205 if ( (int)nx != nx_ ) {
206 d[ 0] = dx_/3.0;
207 d[nx-1] = dx_/3.0;
208 }
209 linear_solve(Mu,dl,d,du,u);
210 }
211
212 void test_inverse_mass(std::ostream &outStream = std::cout) {
213 // State Mass Matrix
214 std::vector<Real> u(nx_,0.0), Mu(nx_,0.0), iMMu(nx_,0.0), diff(nx_,0.0);
215 for (int i = 0; i < nx_; i++) {
216 u[i] = 2.0*(Real)rand()/(Real)RAND_MAX - 1.0;
217 }
218 apply_mass(Mu,u);
219 apply_inverse_mass(iMMu,Mu);
220 axpy(diff,-1.0,iMMu,u);
221 Real error = compute_L2_norm(diff);
222 Real normu = compute_L2_norm(u);
223 outStream << "Test Inverse State Mass Matrix\n";
224 outStream << " ||u - inv(M)Mu|| = " << error << "\n";
225 outStream << " ||u|| = " << normu << "\n";
226 outStream << " Relative Error = " << error/normu << "\n";
227 outStream << "\n";
228 // Control Mass Matrix
229 u.resize(nx_+2,0.0); Mu.resize(nx_+2,0.0); iMMu.resize(nx_+2,0.0); diff.resize(nx_+2,0.0);
230 for (int i = 0; i < nx_+2; i++) {
231 u[i] = 2.0*(Real)rand()/(Real)RAND_MAX - 1.0;
232 }
233 apply_mass(Mu,u);
234 apply_inverse_mass(iMMu,Mu);
235 axpy(diff,-1.0,iMMu,u);
236 error = compute_L2_norm(diff);
237 normu = compute_L2_norm(u);
238 outStream << "Test Inverse Control Mass Matrix\n";
239 outStream << " ||z - inv(M)Mz|| = " << error << "\n";
240 outStream << " ||z|| = " << normu << "\n";
241 outStream << " Relative Error = " << error/normu << "\n";
242 outStream << "\n";
243 }
244
245 /***************************************************************************/
246 /*********************** H1 INNER PRODUCT FUNCTIONS ************************/
247 /***************************************************************************/
248 // Compute H1 inner product
249 Real compute_H1_dot(const std::vector<Real> &x, const std::vector<Real> &y) const {
250 Real ip = 0.0;
251 for (int i=0; i<nx_; i++) {
252 if ( i == 0 ) {
253 ip += cL2_*dx_/6.0*(4.0*x[i] + x[i+1])*y[i]; // Mass term
254 ip += cH1_*(2.0*x[i] - x[i+1])/dx_*y[i]; // Stiffness term
255 }
256 else if ( i == nx_-1 ) {
257 ip += cL2_*dx_/6.0*(x[i-1] + 4.0*x[i])*y[i]; // Mass term
258 ip += cH1_*(2.0*x[i] - x[i-1])/dx_*y[i]; // Stiffness term
259 }
260 else {
261 ip += cL2_*dx_/6.0*(x[i-1] + 4.0*x[i] + x[i+1])*y[i]; // Mass term
262 ip += cH1_*(2.0*x[i] - x[i-1] - x[i+1])/dx_*y[i]; // Stiffness term
263 }
264 }
265 return ip;
266 }
267
268 // compute H1 norm
269 Real compute_H1_norm(const std::vector<Real> &r) const {
270 return std::sqrt(compute_H1_dot(r,r));
271 }
272
273 // Apply H2 Reisz operator
274 void apply_H1(std::vector<Real> &Mu, const std::vector<Real> &u ) const {
275 Mu.resize(nx_,0.0);
276 for (int i=0; i<nx_; i++) {
277 if ( i == 0 ) {
278 Mu[i] = cL2_*dx_/6.0*(4.0*u[i] + u[i+1])
279 + cH1_*(2.0*u[i] - u[i+1])/dx_;
280 }
281 else if ( i == nx_-1 ) {
282 Mu[i] = cL2_*dx_/6.0*(u[i-1] + 4.0*u[i])
283 + cH1_*(2.0*u[i] - u[i-1])/dx_;
284 }
285 else {
286 Mu[i] = cL2_*dx_/6.0*(u[i-1] + 4.0*u[i] + u[i+1])
287 + cH1_*(2.0*u[i] - u[i-1] - u[i+1])/dx_;
288 }
289 }
290 }
291
292 // Apply H1 inverse Reisz operator
293 void apply_inverse_H1(std::vector<Real> &Mu, const std::vector<Real> &u) const {
294 // Build mass matrix
295 std::vector<Real> dl(nx_-1,cL2_*dx_/6.0 - cH1_/dx_);
296 std::vector<Real> d(nx_,2.0*(cL2_*dx_/3.0 + cH1_/dx_));
297 std::vector<Real> du(nx_-1,cL2_*dx_/6.0 - cH1_/dx_);
298 linear_solve(Mu,dl,d,du,u);
299 }
300
301 void test_inverse_H1(std::ostream &outStream = std::cout) {
302 std::vector<Real> u(nx_,0.0), Mu(nx_,0.0), iMMu(nx_,0.0), diff(nx_,0.0);
303 for (int i = 0; i < nx_; i++) {
304 u[i] = 2.0*(Real)rand()/(Real)RAND_MAX - 1.0;
305 }
306 apply_H1(Mu,u);
307 apply_inverse_H1(iMMu,Mu);
308 axpy(diff,-1.0,iMMu,u);
309 Real error = compute_H1_norm(diff);
310 Real normu = compute_H1_norm(u);
311 outStream << "Test Inverse State H1 Matrix\n";
312 outStream << " ||u - inv(M)Mu|| = " << error << "\n";
313 outStream << " ||u|| = " << normu << "\n";
314 outStream << " Relative Error = " << error/normu << "\n";
315 outStream << "\n";
316 }
317
318 /***************************************************************************/
319 /*********************** PDE RESIDUAL AND SOLVE ****************************/
320 /***************************************************************************/
321 // Compute current PDE residual
322 void compute_residual(std::vector<Real> &r, const std::vector<Real> &u,
323 const std::vector<Real> &z) const {
324 r.clear();
325 r.resize(nx_,0.0);
326 for (int i=0; i<nx_; i++) {
327 // Contribution from stiffness term
328 if ( i==0 ) {
329 r[i] = nu_/dx_*(2.0*u[i]-u[i+1]);
330 }
331 else if (i==nx_-1) {
332 r[i] = nu_/dx_*(2.0*u[i]-u[i-1]);
333 }
334 else {
335 r[i] = nu_/dx_*(2.0*u[i]-u[i-1]-u[i+1]);
336 }
337 // Contribution from nonlinear term
338 if (i<nx_-1){
339 r[i] += nl_*u[i+1]*(u[i]+u[i+1])/6.0;
340 }
341 if (i>0) {
342 r[i] -= nl_*u[i-1]*(u[i-1]+u[i])/6.0;
343 }
344 // Contribution from control
345 r[i] -= dx_/6.0*(z[i]+4.0*z[i+1]+z[i+2]);
346 // Contribution from right-hand side
347 r[i] -= dx_*f_;
348 }
349 // Contribution from Dirichlet boundary terms
350 r[0] -= nl_*(u0_*u[ 0]/6.0 + u0_*u0_/6.0) + nu_*u0_/dx_;
351 r[nx_-1] += nl_*(u1_*u[nx_-1]/6.0 + u1_*u1_/6.0) - nu_*u1_/dx_;
352 }
353
354 /***************************************************************************/
355 /*********************** PDE JACOBIAN FUNCTIONS ****************************/
356 /***************************************************************************/
357 // Build PDE Jacobian trigiagonal matrix
358 void compute_pde_jacobian(std::vector<Real> &dl, // Lower diagonal
359 std::vector<Real> &d, // Diagonal
360 std::vector<Real> &du, // Upper diagonal
361 const std::vector<Real> &u) const { // State variable
362 // Get Diagonal and Off-Diagonal Entries of linear PDE Jacobian
363 d.clear();
364 d.resize(nx_,nu_*2.0/dx_);
365 dl.clear();
366 dl.resize(nx_-1,-nu_/dx_);
367 du.clear();
368 du.resize(nx_-1,-nu_/dx_);
369 // Contribution from nonlinearity
370 for (int i=0; i<nx_; i++) {
371 if (i<nx_-1) {
372 dl[i] += nl_*(-2.0*u[i]-u[i+1])/6.0;
373 d[i] += nl_*u[i+1]/6.0;
374 }
375 if (i>0) {
376 d[i] -= nl_*u[i-1]/6.0;
377 du[i-1] += nl_*(u[i-1]+2.0*u[i])/6.0;
378 }
379 }
380 // Contribution from Dirichlet boundary conditions
381 d[ 0] -= nl_*u0_/6.0;
382 d[nx_-1] += nl_*u1_/6.0;
383 }
384
385 // Apply PDE Jacobian to a vector
386 void apply_pde_jacobian(std::vector<Real> &jv,
387 const std::vector<Real> &v,
388 const std::vector<Real> &u,
389 const std::vector<Real> &z) const {
390 // Fill jv
391 for (int i = 0; i < nx_; i++) {
392 jv[i] = nu_/dx_*2.0*v[i];
393 if ( i > 0 ) {
394 jv[i] += -nu_/dx_*v[i-1]-nl_*(u[i-1]/6.0*v[i]+(u[i]+2.0*u[i-1])/6.0*v[i-1]);
395 }
396 if ( i < nx_-1 ) {
397 jv[i] += -nu_/dx_*v[i+1]+nl_*(u[i+1]/6.0*v[i]+(u[i]+2.0*u[i+1])/6.0*v[i+1]);
398 }
399 }
400 jv[ 0] -= nl_*u0_/6.0*v[0];
401 jv[nx_-1] += nl_*u1_/6.0*v[nx_-1];
402 }
403
404 // Apply inverse PDE Jacobian to a vector
405 void apply_inverse_pde_jacobian(std::vector<Real> &ijv,
406 const std::vector<Real> &v,
407 const std::vector<Real> &u,
408 const std::vector<Real> &z) const {
409 // Get PDE Jacobian
410 std::vector<Real> d(nx_,0.0);
411 std::vector<Real> dl(nx_-1,0.0);
412 std::vector<Real> du(nx_-1,0.0);
413 compute_pde_jacobian(dl,d,du,u);
414 // Solve solve state sensitivity system at current time step
415 linear_solve(ijv,dl,d,du,v);
416 }
417
418 // Apply adjoint PDE Jacobian to a vector
419 void apply_adjoint_pde_jacobian(std::vector<Real> &ajv,
420 const std::vector<Real> &v,
421 const std::vector<Real> &u,
422 const std::vector<Real> &z) const {
423 // Fill jvp
424 for (int i = 0; i < nx_; i++) {
425 ajv[i] = nu_/dx_*2.0*v[i];
426 if ( i > 0 ) {
427 ajv[i] += -nu_/dx_*v[i-1]-nl_*(u[i-1]/6.0*v[i]
428 -(u[i-1]+2.0*u[i])/6.0*v[i-1]);
429 }
430 if ( i < nx_-1 ) {
431 ajv[i] += -nu_/dx_*v[i+1]+nl_*(u[i+1]/6.0*v[i]
432 -(u[i+1]+2.0*u[i])/6.0*v[i+1]);
433 }
434 }
435 ajv[ 0] -= nl_*u0_/6.0*v[0];
436 ajv[nx_-1] += nl_*u1_/6.0*v[nx_-1];
437 }
438
439 // Apply inverse adjoint PDE Jacobian to a vector
440 void apply_inverse_adjoint_pde_jacobian(std::vector<Real> &iajv,
441 const std::vector<Real> &v,
442 const std::vector<Real> &u,
443 const std::vector<Real> &z) const {
444 // Get PDE Jacobian
445 std::vector<Real> d(nx_,0.0);
446 std::vector<Real> du(nx_-1,0.0);
447 std::vector<Real> dl(nx_-1,0.0);
448 compute_pde_jacobian(dl,d,du,u);
449 // Solve solve adjoint system at current time step
450 linear_solve(iajv,dl,d,du,v,true);
451 }
452
453 /***************************************************************************/
454 /*********************** CONTROL JACOBIAN FUNCTIONS ************************/
455 /***************************************************************************/
456 // Apply control Jacobian to a vector
457 void apply_control_jacobian(std::vector<Real> &jv,
458 const std::vector<Real> &v,
459 const std::vector<Real> &u,
460 const std::vector<Real> &z) const {
461 for (int i=0; i<nx_; i++) {
462 // Contribution from control
463 jv[i] = -dx_/6.0*(v[i]+4.0*v[i+1]+v[i+2]);
464 }
465 }
466
467 // Apply adjoint control Jacobian to a vector
468 void apply_adjoint_control_jacobian(std::vector<Real> &jv,
469 const std::vector<Real> &v,
470 const std::vector<Real> &u,
471 const std::vector<Real> &z) const {
472 for (int i=0; i<nx_+2; i++) {
473 if ( i == 0 ) {
474 jv[i] = -dx_/6.0*v[i];
475 }
476 else if ( i == 1 ) {
477 jv[i] = -dx_/6.0*(4.0*v[i-1]+v[i]);
478 }
479 else if ( i == nx_ ) {
480 jv[i] = -dx_/6.0*(4.0*v[i-1]+v[i-2]);
481 }
482 else if ( i == nx_+1 ) {
483 jv[i] = -dx_/6.0*v[i-2];
484 }
485 else {
486 jv[i] = -dx_/6.0*(v[i-2]+4.0*v[i-1]+v[i]);
487 }
488 }
489 }
490
491 /***************************************************************************/
492 /*********************** AJDOINT HESSIANS **********************************/
493 /***************************************************************************/
494 void apply_adjoint_pde_hessian(std::vector<Real> &ahwv,
495 const std::vector<Real> &w,
496 const std::vector<Real> &v,
497 const std::vector<Real> &u,
498 const std::vector<Real> &z) const {
499 for (int i=0; i<nx_; i++) {
500 // Contribution from nonlinear term
501 ahwv[i] = 0.0;
502 if (i<nx_-1){
503 ahwv[i] += (w[i]*v[i+1] - w[i+1]*(2.0*v[i]+v[i+1]))/6.0;
504 }
505 if (i>0) {
506 ahwv[i] += (w[i-1]*(v[i-1]+2.0*v[i]) - w[i]*v[i-1])/6.0;
507 }
508 }
509 //ahwv.assign(u.size(),0.0);
510 }
511 void apply_adjoint_pde_control_hessian(std::vector<Real> &ahwv,
512 const std::vector<Real> &w,
513 const std::vector<Real> &v,
514 const std::vector<Real> &u,
515 const std::vector<Real> &z) {
516 ahwv.assign(u.size(),0.0);
517 }
518 void apply_adjoint_control_pde_hessian(std::vector<Real> &ahwv,
519 const std::vector<Real> &w,
520 const std::vector<Real> &v,
521 const std::vector<Real> &u,
522 const std::vector<Real> &z) {
523 ahwv.assign(z.size(),0.0);
524 }
525 void apply_adjoint_control_hessian(std::vector<Real> &ahwv,
526 const std::vector<Real> &w,
527 const std::vector<Real> &v,
528 const std::vector<Real> &u,
529 const std::vector<Real> &z) {
530 ahwv.assign(z.size(),0.0);
531 }
532};
533
534template<class Real>
535class L2VectorPrimal : public ROL::Vector<Real> {
536private:
537 ROL::Ptr<std::vector<Real> > vec_;
538 ROL::Ptr<BurgersFEM<Real> > fem_;
539
540 mutable ROL::Ptr<L2VectorDual<Real> > dual_vec_;
541
542public:
543 L2VectorPrimal(const ROL::Ptr<std::vector<Real> > & vec,
544 const ROL::Ptr<BurgersFEM<Real> > &fem)
545 : vec_(vec), fem_(fem), dual_vec_(ROL::nullPtr) {}
546
547 void set( const ROL::Vector<Real> &x ) {
548 const L2VectorPrimal &ex = dynamic_cast<const L2VectorPrimal&>(x);
549 const std::vector<Real>& xval = *ex.getVector();
550 std::copy(xval.begin(),xval.end(),vec_->begin());
551 }
552
553 void plus( const ROL::Vector<Real> &x ) {
554 const L2VectorPrimal &ex = dynamic_cast<const L2VectorPrimal&>(x);
555 const std::vector<Real>& xval = *ex.getVector();
556 unsigned dimension = vec_->size();
557 for (unsigned i=0; i<dimension; i++) {
558 (*vec_)[i] += xval[i];
559 }
560 }
561
562 void scale( const Real alpha ) {
563 unsigned dimension = vec_->size();
564 for (unsigned i=0; i<dimension; i++) {
565 (*vec_)[i] *= alpha;
566 }
567 }
568
569 Real dot( const ROL::Vector<Real> &x ) const {
570 const L2VectorPrimal & ex = dynamic_cast<const L2VectorPrimal&>(x);
571 const std::vector<Real>& xval = *ex.getVector();
572 return fem_->compute_L2_dot(xval,*vec_);
573 }
574
575 Real norm() const {
576 Real val = 0;
577 val = std::sqrt( dot(*this) );
578 return val;
579 }
580
581 ROL::Ptr<ROL::Vector<Real> > clone() const {
582 return ROL::makePtr<L2VectorPrimal>( ROL::makePtr<std::vector<Real>>(vec_->size(),0.0),fem_);
583 }
584
585 ROL::Ptr<const std::vector<Real> > getVector() const {
586 return vec_;
587 }
588
589 ROL::Ptr<std::vector<Real> > getVector() {
590 return vec_;
591 }
592
593 ROL::Ptr<ROL::Vector<Real> > basis( const int i ) const {
594 ROL::Ptr<L2VectorPrimal> e
595 = ROL::makePtr<L2VectorPrimal>( ROL::makePtr<std::vector<Real>>(vec_->size(),0.0),fem_);
596 (*e->getVector())[i] = 1.0;
597 return e;
598 }
599
600 int dimension() const {
601 return vec_->size();
602 }
603
604 const ROL::Vector<Real>& dual() const {
605 dual_vec_ = ROL::makePtr<L2VectorDual<Real>>(
606 ROL::makePtr<std::vector<Real>>(*vec_),fem_);
607
608 fem_->apply_mass(*(ROL::constPtrCast<std::vector<Real> >(dual_vec_->getVector())),*vec_);
609 return *dual_vec_;
610 }
611
612 Real apply(const ROL::Vector<Real> &x) const {
613 const L2VectorDual<Real> &ex = dynamic_cast<const L2VectorDual<Real>&>(x);
614 const std::vector<Real>& xval = *ex.getVector();
615 return std::inner_product(vec_->begin(), vec_->end(), xval.begin(), Real(0));
616 }
617
618};
619
620template<class Real>
621class L2VectorDual : public ROL::Vector<Real> {
622private:
623 ROL::Ptr<std::vector<Real> > vec_;
624 ROL::Ptr<BurgersFEM<Real> > fem_;
625
626 mutable ROL::Ptr<L2VectorPrimal<Real> > dual_vec_;
627
628public:
629 L2VectorDual(const ROL::Ptr<std::vector<Real> > & vec,
630 const ROL::Ptr<BurgersFEM<Real> > &fem)
631 : vec_(vec), fem_(fem), dual_vec_(ROL::nullPtr) {}
632
633 void set( const ROL::Vector<Real> &x ) {
634 const L2VectorDual &ex = dynamic_cast<const L2VectorDual&>(x);
635 const std::vector<Real>& xval = *ex.getVector();
636 std::copy(xval.begin(),xval.end(),vec_->begin());
637 }
638
639 void plus( const ROL::Vector<Real> &x ) {
640 const L2VectorDual &ex = dynamic_cast<const L2VectorDual&>(x);
641 const std::vector<Real>& xval = *ex.getVector();
642 unsigned dimension = vec_->size();
643 for (unsigned i=0; i<dimension; i++) {
644 (*vec_)[i] += xval[i];
645 }
646 }
647
648 void scale( const Real alpha ) {
649 unsigned dimension = vec_->size();
650 for (unsigned i=0; i<dimension; i++) {
651 (*vec_)[i] *= alpha;
652 }
653 }
654
655 Real dot( const ROL::Vector<Real> &x ) const {
656 const L2VectorDual & ex = dynamic_cast<const L2VectorDual&>(x);
657 const std::vector<Real>& xval = *ex.getVector();
658 unsigned dimension = vec_->size();
659 std::vector<Real> Mx(dimension,0.0);
660 fem_->apply_inverse_mass(Mx,xval);
661 Real val = 0.0;
662 for (unsigned i = 0; i < dimension; i++) {
663 val += Mx[i]*(*vec_)[i];
664 }
665 return val;
666 }
667
668 Real norm() const {
669 Real val = 0;
670 val = std::sqrt( dot(*this) );
671 return val;
672 }
673
674 ROL::Ptr<ROL::Vector<Real> > clone() const {
675 return ROL::makePtr<L2VectorDual>( ROL::makePtr<std::vector<Real>>(vec_->size(),0.0),fem_);
676 }
677
678 ROL::Ptr<const std::vector<Real> > getVector() const {
679 return vec_;
680 }
681
682 ROL::Ptr<std::vector<Real> > getVector() {
683 return vec_;
684 }
685
686 ROL::Ptr<ROL::Vector<Real> > basis( const int i ) const {
687 ROL::Ptr<L2VectorDual> e
688 = ROL::makePtr<L2VectorDual>( ROL::makePtr<std::vector<Real>>(vec_->size(),0.0),fem_);
689 (*e->getVector())[i] = 1.0;
690 return e;
691 }
692
693 int dimension() const {
694 return vec_->size();
695 }
696
697 const ROL::Vector<Real>& dual() const {
698 dual_vec_ = ROL::makePtr<L2VectorPrimal<Real>>(
699 ROL::makePtr<std::vector<Real>>(*vec_),fem_);
700
701 fem_->apply_inverse_mass(*(ROL::constPtrCast<std::vector<Real> >(dual_vec_->getVector())),*vec_);
702 return *dual_vec_;
703 }
704
705 Real apply(const ROL::Vector<Real> &x) const {
706 const L2VectorPrimal<Real> &ex = dynamic_cast<const L2VectorPrimal<Real>&>(x);
707 const std::vector<Real>& xval = *ex.getVector();
708 return std::inner_product(vec_->begin(), vec_->end(), xval.begin(), Real(0));
709 }
710
711};
712
713template<class Real>
714class H1VectorPrimal : public ROL::Vector<Real> {
715private:
716 ROL::Ptr<std::vector<Real> > vec_;
717 ROL::Ptr<BurgersFEM<Real> > fem_;
718
719 mutable ROL::Ptr<H1VectorDual<Real> > dual_vec_;
720
721public:
722 H1VectorPrimal(const ROL::Ptr<std::vector<Real> > & vec,
723 const ROL::Ptr<BurgersFEM<Real> > &fem)
724 : vec_(vec), fem_(fem), dual_vec_(ROL::nullPtr) {}
725
726 void set( const ROL::Vector<Real> &x ) {
727 const H1VectorPrimal &ex = dynamic_cast<const H1VectorPrimal&>(x);
728 const std::vector<Real>& xval = *ex.getVector();
729 std::copy(xval.begin(),xval.end(),vec_->begin());
730 }
731
732 void plus( const ROL::Vector<Real> &x ) {
733 const H1VectorPrimal &ex = dynamic_cast<const H1VectorPrimal&>(x);
734 const std::vector<Real>& xval = *ex.getVector();
735 unsigned dimension = vec_->size();
736 for (unsigned i=0; i<dimension; i++) {
737 (*vec_)[i] += xval[i];
738 }
739 }
740
741 void scale( const Real alpha ) {
742 unsigned dimension = vec_->size();
743 for (unsigned i=0; i<dimension; i++) {
744 (*vec_)[i] *= alpha;
745 }
746 }
747
748 Real dot( const ROL::Vector<Real> &x ) const {
749 const H1VectorPrimal & ex = dynamic_cast<const H1VectorPrimal&>(x);
750 const std::vector<Real>& xval = *ex.getVector();
751 return fem_->compute_H1_dot(xval,*vec_);
752 }
753
754 Real norm() const {
755 Real val = 0;
756 val = std::sqrt( dot(*this) );
757 return val;
758 }
759
760 ROL::Ptr<ROL::Vector<Real> > clone() const {
761 return ROL::makePtr<H1VectorPrimal>( ROL::makePtr<std::vector<Real>>(vec_->size(),0.0),fem_);
762 }
763
764 ROL::Ptr<const std::vector<Real> > getVector() const {
765 return vec_;
766 }
767
768 ROL::Ptr<std::vector<Real> > getVector() {
769 return vec_;
770 }
771
772 ROL::Ptr<ROL::Vector<Real> > basis( const int i ) const {
773 ROL::Ptr<H1VectorPrimal> e
774 = ROL::makePtr<H1VectorPrimal>( ROL::makePtr<std::vector<Real>>(vec_->size(),0.0),fem_);
775 (*e->getVector())[i] = 1.0;
776 return e;
777 }
778
779 int dimension() const {
780 return vec_->size();
781 }
782
783 const ROL::Vector<Real>& dual() const {
784 dual_vec_ = ROL::makePtr<H1VectorDual<Real>>(
785 ROL::makePtr<std::vector<Real>>(*vec_),fem_);
786
787 fem_->apply_H1(*(ROL::constPtrCast<std::vector<Real> >(dual_vec_->getVector())),*vec_);
788 return *dual_vec_;
789 }
790
791 Real apply(const ROL::Vector<Real> &x) const {
792 const H1VectorDual<Real> &ex = dynamic_cast<const H1VectorDual<Real>&>(x);
793 const std::vector<Real>& xval = *ex.getVector();
794 return std::inner_product(vec_->begin(), vec_->end(), xval.begin(), Real(0));
795 }
796
797};
798
799template<class Real>
800class H1VectorDual : public ROL::Vector<Real> {
801private:
802 ROL::Ptr<std::vector<Real> > vec_;
803 ROL::Ptr<BurgersFEM<Real> > fem_;
804
805 mutable ROL::Ptr<H1VectorPrimal<Real> > dual_vec_;
806
807public:
808 H1VectorDual(const ROL::Ptr<std::vector<Real> > & vec,
809 const ROL::Ptr<BurgersFEM<Real> > &fem)
810 : vec_(vec), fem_(fem), dual_vec_(ROL::nullPtr) {}
811
812 void set( const ROL::Vector<Real> &x ) {
813 const H1VectorDual &ex = dynamic_cast<const H1VectorDual&>(x);
814 const std::vector<Real>& xval = *ex.getVector();
815 std::copy(xval.begin(),xval.end(),vec_->begin());
816 }
817
818 void plus( const ROL::Vector<Real> &x ) {
819 const H1VectorDual &ex = dynamic_cast<const H1VectorDual&>(x);
820 const std::vector<Real>& xval = *ex.getVector();
821 unsigned dimension = vec_->size();
822 for (unsigned i=0; i<dimension; i++) {
823 (*vec_)[i] += xval[i];
824 }
825 }
826
827 void scale( const Real alpha ) {
828 unsigned dimension = vec_->size();
829 for (unsigned i=0; i<dimension; i++) {
830 (*vec_)[i] *= alpha;
831 }
832 }
833
834 Real dot( const ROL::Vector<Real> &x ) const {
835 const H1VectorDual & ex = dynamic_cast<const H1VectorDual&>(x);
836 const std::vector<Real>& xval = *ex.getVector();
837 unsigned dimension = vec_->size();
838 std::vector<Real> Mx(dimension,0.0);
839 fem_->apply_inverse_H1(Mx,xval);
840 Real val = 0.0;
841 for (unsigned i = 0; i < dimension; i++) {
842 val += Mx[i]*(*vec_)[i];
843 }
844 return val;
845 }
846
847 Real norm() const {
848 Real val = 0;
849 val = std::sqrt( dot(*this) );
850 return val;
851 }
852
853 ROL::Ptr<ROL::Vector<Real> > clone() const {
854 return ROL::makePtr<H1VectorDual>( ROL::makePtr<std::vector<Real>>(vec_->size(),0.0),fem_);
855 }
856
857 ROL::Ptr<const std::vector<Real> > getVector() const {
858 return vec_;
859 }
860
861 ROL::Ptr<std::vector<Real> > getVector() {
862 return vec_;
863 }
864
865 ROL::Ptr<ROL::Vector<Real> > basis( const int i ) const {
866 ROL::Ptr<H1VectorDual> e
867 = ROL::makePtr<H1VectorDual>( ROL::makePtr<std::vector<Real>>(vec_->size(),0.0),fem_);
868 (*e->getVector())[i] = 1.0;
869 return e;
870 }
871
872 int dimension() const {
873 return vec_->size();
874 }
875
876 const ROL::Vector<Real>& dual() const {
877 dual_vec_ = ROL::makePtr<H1VectorPrimal<Real>>(
878 ROL::makePtr<std::vector<Real>>(*vec_),fem_);
879
880 fem_->apply_inverse_H1(*(ROL::constPtrCast<std::vector<Real> >(dual_vec_->getVector())),*vec_);
881 return *dual_vec_;
882 }
883
884 Real apply(const ROL::Vector<Real> &x) const {
885 const H1VectorPrimal<Real> &ex = dynamic_cast<const H1VectorPrimal<Real>&>(x);
886 const std::vector<Real>& xval = *ex.getVector();
887 return std::inner_product(vec_->begin(), vec_->end(), xval.begin(), Real(0));
888 }
889};
890
891
892template<class Real>
894private:
895
898
901
904
905 typedef typename std::vector<Real>::size_type uint;
906
907 ROL::Ptr<BurgersFEM<Real> > fem_;
908 bool useHessian_;
909
910public:
911 Constraint_BurgersControl(ROL::Ptr<BurgersFEM<Real> > &fem, bool useHessian = true)
912 : fem_(fem), useHessian_(useHessian) {}
913
915 const ROL::Vector<Real> &z, Real &tol) {
916 ROL::Ptr<std::vector<Real> > cp =
917 dynamic_cast<PrimalConstraintVector&>(c).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_->compute_residual(*cp,*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 PrimalStateVector&>(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
941 const std::vector<Real> param
943 fem_->set_problem_data(param[0],param[1],param[2],param[3]);
944
945 fem_->apply_pde_jacobian(*jvp,*vp,*up,*zp);
946 }
947
949 const ROL::Vector<Real> &z, Real &tol) {
950 ROL::Ptr<std::vector<Real> > jvp =
951 dynamic_cast<PrimalConstraintVector&>(jv).getVector();
952 ROL::Ptr<const std::vector<Real> > vp =
953 dynamic_cast<const PrimalControlVector&>(v).getVector();
954 ROL::Ptr<const std::vector<Real> > up =
955 dynamic_cast<const PrimalStateVector&>(u).getVector();
956 ROL::Ptr<const std::vector<Real> > zp =
957 dynamic_cast<const PrimalControlVector&>(z).getVector();
958
959 const std::vector<Real> param
961 fem_->set_problem_data(param[0],param[1],param[2],param[3]);
962
963 fem_->apply_control_jacobian(*jvp,*vp,*up,*zp);
964 }
965
967 const ROL::Vector<Real> &z, Real &tol) {
968 ROL::Ptr<std::vector<Real> > ijvp =
969 dynamic_cast<PrimalStateVector&>(ijv).getVector();
970 ROL::Ptr<const std::vector<Real> > vp =
971 dynamic_cast<const PrimalConstraintVector&>(v).getVector();
972 ROL::Ptr<const std::vector<Real> > up =
973 dynamic_cast<const PrimalStateVector&>(u).getVector();
974 ROL::Ptr<const std::vector<Real> > zp =
975 dynamic_cast<const PrimalControlVector&>(z).getVector();
976
977 const std::vector<Real> param
979 fem_->set_problem_data(param[0],param[1],param[2],param[3]);
980
981 fem_->apply_inverse_pde_jacobian(*ijvp,*vp,*up,*zp);
982 }
983
985 const ROL::Vector<Real> &z, Real &tol) {
986 ROL::Ptr<std::vector<Real> > jvp =
987 dynamic_cast<DualStateVector&>(ajv).getVector();
988 ROL::Ptr<const std::vector<Real> > vp =
989 dynamic_cast<const DualConstraintVector&>(v).getVector();
990 ROL::Ptr<const std::vector<Real> > up =
991 dynamic_cast<const PrimalStateVector&>(u).getVector();
992 ROL::Ptr<const std::vector<Real> > zp =
993 dynamic_cast<const PrimalControlVector&>(z).getVector();
994
995 const std::vector<Real> param
997 fem_->set_problem_data(param[0],param[1],param[2],param[3]);
998
999 fem_->apply_adjoint_pde_jacobian(*jvp,*vp,*up,*zp);
1000 }
1001
1003 const ROL::Vector<Real> &z, Real &tol) {
1004 ROL::Ptr<std::vector<Real> > jvp =
1005 dynamic_cast<DualControlVector&>(jv).getVector();
1006 ROL::Ptr<const std::vector<Real> > vp =
1007 dynamic_cast<const DualConstraintVector&>(v).getVector();
1008 ROL::Ptr<const std::vector<Real> > up =
1009 dynamic_cast<const PrimalStateVector&>(u).getVector();
1010 ROL::Ptr<const std::vector<Real> > zp =
1011 dynamic_cast<const PrimalControlVector&>(z).getVector();
1012
1013 const std::vector<Real> param
1015 fem_->set_problem_data(param[0],param[1],param[2],param[3]);
1016
1017 fem_->apply_adjoint_control_jacobian(*jvp,*vp,*up,*zp);
1018 }
1019
1021 const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
1022 ROL::Ptr<std::vector<Real> > iajvp =
1023 dynamic_cast<DualConstraintVector&>(iajv).getVector();
1024 ROL::Ptr<const std::vector<Real> > vp =
1025 dynamic_cast<const DualStateVector&>(v).getVector();
1026 ROL::Ptr<const std::vector<Real> > up =
1027 dynamic_cast<const PrimalStateVector&>(u).getVector();
1028 ROL::Ptr<const std::vector<Real> > zp =
1029 dynamic_cast<const PrimalControlVector&>(z).getVector();
1030
1031 const std::vector<Real> param
1033 fem_->set_problem_data(param[0],param[1],param[2],param[3]);
1034
1035 fem_->apply_inverse_adjoint_pde_jacobian(*iajvp,*vp,*up,*zp);
1036 }
1037
1039 const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
1040 if ( useHessian_ ) {
1041 ROL::Ptr<std::vector<Real> > ahwvp =
1042 dynamic_cast<DualStateVector&>(ahwv).getVector();
1043 ROL::Ptr<const std::vector<Real> > wp =
1044 dynamic_cast<const DualConstraintVector&>(w).getVector();
1045 ROL::Ptr<const std::vector<Real> > vp =
1046 dynamic_cast<const PrimalStateVector&>(v).getVector();
1047 ROL::Ptr<const std::vector<Real> > up =
1048 dynamic_cast<const PrimalStateVector&>(u).getVector();
1049 ROL::Ptr<const std::vector<Real> > zp =
1050 dynamic_cast<const PrimalControlVector&>(z).getVector();
1051
1052 const std::vector<Real> param
1054 fem_->set_problem_data(param[0],param[1],param[2],param[3]);
1055
1056 fem_->apply_adjoint_pde_hessian(*ahwvp,*wp,*vp,*up,*zp);
1057 }
1058 else {
1059 ahwv.zero();
1060 }
1061 }
1062
1064 const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
1065 if ( useHessian_ ) {
1066 ROL::Ptr<std::vector<Real> > ahwvp =
1067 dynamic_cast<DualControlVector&>(ahwv).getVector();
1068 ROL::Ptr<const std::vector<Real> > wp =
1069 dynamic_cast<const DualConstraintVector&>(w).getVector();
1070 ROL::Ptr<const std::vector<Real> > vp =
1071 dynamic_cast<const PrimalStateVector&>(v).getVector();
1072 ROL::Ptr<const std::vector<Real> > up =
1073 dynamic_cast<const PrimalStateVector&>(u).getVector();
1074 ROL::Ptr<const std::vector<Real> > zp =
1075 dynamic_cast<const PrimalControlVector&>(z).getVector();
1076
1077 const std::vector<Real> param
1079 fem_->set_problem_data(param[0],param[1],param[2],param[3]);
1080
1081 fem_->apply_adjoint_control_pde_hessian(*ahwvp,*wp,*vp,*up,*zp);
1082 }
1083 else {
1084 ahwv.zero();
1085 }
1086 }
1088 const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
1089 if ( useHessian_ ) {
1090 ROL::Ptr<std::vector<Real> > ahwvp =
1091 dynamic_cast<DualStateVector&>(ahwv).getVector();
1092 ROL::Ptr<const std::vector<Real> > wp =
1093 dynamic_cast<const DualConstraintVector&>(w).getVector();
1094 ROL::Ptr<const std::vector<Real> > vp =
1095 dynamic_cast<const PrimalControlVector&>(v).getVector();
1096 ROL::Ptr<const std::vector<Real> > up =
1097 dynamic_cast<const PrimalStateVector&>(u).getVector();
1098 ROL::Ptr<const std::vector<Real> > zp =
1099 dynamic_cast<const PrimalControlVector&>(z).getVector();
1100
1101 const std::vector<Real> param
1103 fem_->set_problem_data(param[0],param[1],param[2],param[3]);
1104
1105 fem_->apply_adjoint_pde_control_hessian(*ahwvp,*wp,*vp,*up,*zp);
1106 }
1107 else {
1108 ahwv.zero();
1109 }
1110 }
1112 const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
1113 if ( useHessian_ ) {
1114 ROL::Ptr<std::vector<Real> > ahwvp =
1115 dynamic_cast<DualControlVector&>(ahwv).getVector();
1116 ROL::Ptr<const std::vector<Real> > wp =
1117 dynamic_cast<const DualConstraintVector&>(w).getVector();
1118 ROL::Ptr<const std::vector<Real> > vp =
1119 dynamic_cast<const PrimalControlVector&>(v).getVector();
1120 ROL::Ptr<const std::vector<Real> > up =
1121 dynamic_cast<const PrimalStateVector&>(u).getVector();
1122 ROL::Ptr<const std::vector<Real> > zp =
1123 dynamic_cast<const PrimalControlVector&>(z).getVector();
1124
1125 const std::vector<Real> param
1127 fem_->set_problem_data(param[0],param[1],param[2],param[3]);
1128
1129 fem_->apply_adjoint_control_hessian(*ahwvp,*wp,*vp,*up,*zp);
1130 }
1131 else {
1132 ahwv.zero();
1133 }
1134 }
1135};
1136
1137template<class Real>
1139private:
1140
1143
1146
1147 typedef typename std::vector<Real>::size_type uint;
1148
1149 ROL::Ptr<BurgersFEM<Real> > fem_;
1150
1151 Real x_;
1152 std::vector<int> indices_;
1153
1154public:
1156 Real x = 0.0) : fem_(fem), x_(x) {
1157 for (int i = 1; i < fem_->num_dof()+1; i++) {
1158 if ( (Real)i*(fem_->mesh_spacing()) >= x_ ) {
1159 indices_.push_back(i-1);
1160 }
1161 }
1162 }
1163
1164 using ROL::Objective_SimOpt<Real>::value;
1165
1166 Real value( const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
1167 ROL::Ptr<const std::vector<Real> > up =
1168 dynamic_cast<const PrimalStateVector&>(u).getVector();
1169
1170// const std::vector<Real> param
1171// = ROL::Objective_SimOpt<Real>::getParameter();
1172// fem_->set_problem_data(param[0],param[1],param[2],param[3]);
1173// Real nu = fem_->get_viscosity();
1174//
1175// return 0.5*nu*fem_->compute_H1_dot(*up,*up);
1176
1177 Real val = 0.5*((((Real)indices_[0]+1.)*(fem_->mesh_spacing())-x_)
1178 *(x_+(2.-((Real)indices_[0]+1.))*(fem_->mesh_spacing()))/(fem_->mesh_spacing())
1179 +(fem_->mesh_spacing())) * (*up)[indices_[0]];
1180 for (uint i = 1; i < indices_.size(); i++) {
1181 val += (fem_->mesh_spacing())*(*up)[indices_[i]];
1182 }
1183 return -val;
1184 }
1185
1186 void gradient_1( ROL::Vector<Real> &g, const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
1187 ROL::Ptr<std::vector<Real> > gp =
1188 dynamic_cast<DualStateVector&>(g).getVector();
1189 ROL::Ptr<const std::vector<Real> > up =
1190 dynamic_cast<const PrimalStateVector&>(u).getVector();
1191
1192// const std::vector<Real> param
1193// = ROL::Objective_SimOpt<Real>::getParameter();
1194// fem_->set_problem_data(param[0],param[1],param[2],param[3]);
1195// Real nu = fem_->get_viscosity();
1196//
1197// fem_->apply_H1(*gp,*up);
1198// g.scale(nu);
1199
1200 g.zero();
1201 (*gp)[indices_[0]] = -0.5*((((Real)indices_[0]+1.)*(fem_->mesh_spacing())-x_)
1202 *(x_+(2.-((Real)indices_[0]+1.))*(fem_->mesh_spacing()))/(fem_->mesh_spacing())
1203 +(fem_->mesh_spacing()));
1204
1205
1206 for (uint i = 1; i < indices_.size(); i++) {
1207 (*gp)[indices_[i]] = -(fem_->mesh_spacing());
1208 }
1209 }
1210
1211 void gradient_2( ROL::Vector<Real> &g, const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
1212 g.zero();
1213 }
1214
1216 const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
1217// ROL::Ptr<std::vector<Real> > hvp =
1218// ROL::constPtrCast<std::vector<Real> >((dynamic_cast<DualStateVector&>(hv)).getVector());
1219// ROL::Ptr<const std::vector<Real> > vp =
1220// (dynamic_cast<PrimalStateVector>(const_cast<ROL::Vector<Real> &&>(v))).getVector();
1221//
1222// const std::vector<Real> param
1223// = ROL::Objective_SimOpt<Real>::getParameter();
1224// fem_->set_problem_data(param[0],param[1],param[2],param[3]);
1225// Real nu = fem_->get_viscosity();
1226//
1227// fem_->apply_H1(*hvp,*vp);
1228// hv.scale(nu);
1229
1230 hv.zero();
1231 }
1232
1234 const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
1235 hv.zero();
1236 }
1237
1239 const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
1240 hv.zero();
1241 }
1242
1244 const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
1245 hv.zero();
1246 }
1247};
1248
1249template<class Real>
1250class L2BoundConstraint : public ROL::BoundConstraint<Real> {
1251private:
1252 int dim_;
1253 std::vector<Real> x_lo_;
1254 std::vector<Real> x_up_;
1255 Real min_diff_;
1256 Real scale_;
1257 ROL::Ptr<BurgersFEM<Real> > fem_;
1258 ROL::Ptr<ROL::Vector<Real> > l_;
1259 ROL::Ptr<ROL::Vector<Real> > u_;
1260
1261 void cast_vector(ROL::Ptr<std::vector<Real> > &xvec,
1262 ROL::Vector<Real> &x) const {
1263 try {
1264 xvec = dynamic_cast<L2VectorPrimal<Real>&>(x).getVector();
1265 }
1266 catch (std::exception &e) {
1267 xvec = dynamic_cast<L2VectorDual<Real>&>(x).getVector();
1268 }
1269 }
1270
1271 void cast_const_vector(ROL::Ptr<const std::vector<Real> > &xvec,
1272 const ROL::Vector<Real> &x) const {
1273 try {
1274 xvec = dynamic_cast<const L2VectorPrimal<Real>&>(x).getVector();
1275 }
1276 catch (std::exception &e) {
1277 xvec = dynamic_cast<const L2VectorDual<Real>&>(x).getVector();
1278 }
1279 }
1280
1281 void axpy(std::vector<Real> &out, const Real a,
1282 const std::vector<Real> &x, const std::vector<Real> &y) const{
1283 out.resize(dim_,0.0);
1284 for (unsigned i = 0; i < dim_; i++) {
1285 out[i] = a*x[i] + y[i];
1286 }
1287 }
1288
1289 void projection(std::vector<Real> &x) {
1290 for ( int i = 0; i < dim_; i++ ) {
1291 x[i] = std::max(x_lo_[i],std::min(x_up_[i],x[i]));
1292 }
1293 }
1294
1295public:
1296 L2BoundConstraint(std::vector<Real> &l, std::vector<Real> &u,
1297 const ROL::Ptr<BurgersFEM<Real> > &fem, Real scale = 1.0)
1298 : x_lo_(l), x_up_(u), scale_(scale), fem_(fem) {
1299 dim_ = x_lo_.size();
1300 for ( int i = 0; i < dim_; i++ ) {
1301 if ( i == 0 ) {
1302 min_diff_ = x_up_[i] - x_lo_[i];
1303 }
1304 else {
1305 min_diff_ = ( (min_diff_ < (x_up_[i] - x_lo_[i])) ? min_diff_ : (x_up_[i] - x_lo_[i]) );
1306 }
1307 }
1308 min_diff_ *= 0.5;
1309 l_ = ROL::makePtr<L2VectorPrimal<Real>>(
1310 ROL::makePtr<std::vector<Real>>(l), fem);
1311 u_ = ROL::makePtr<L2VectorPrimal<Real>>(
1312 ROL::makePtr<std::vector<Real>>(u), fem);
1313 }
1314
1316 ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1317 bool val = true;
1318 int cnt = 1;
1319 for ( int i = 0; i < dim_; i++ ) {
1320 if ( (*ex)[i] >= x_lo_[i] && (*ex)[i] <= x_up_[i] ) { cnt *= 1; }
1321 else { cnt *= 0; }
1322 }
1323 if ( cnt == 0 ) { val = false; }
1324 return val;
1325 }
1326
1328 ROL::Ptr<std::vector<Real> > ex; cast_vector(ex,x);
1329 projection(*ex);
1330 }
1331
1332 void pruneLowerActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &x, Real eps = Real(0)) {
1333 ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1334 ROL::Ptr<std::vector<Real> > ev; cast_vector(ev,v);
1335 Real epsn = std::min(scale_*eps,min_diff_);
1336 for ( int i = 0; i < dim_; i++ ) {
1337 if ( ((*ex)[i] <= x_lo_[i]+epsn) ) {
1338 (*ev)[i] = 0.0;
1339 }
1340 }
1341 }
1342
1343 void pruneUpperActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &x, Real eps = Real(0)) {
1344 ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1345 ROL::Ptr<std::vector<Real> > ev; cast_vector(ev,v);
1346 Real epsn = std::min(scale_*eps,min_diff_);
1347 for ( int i = 0; i < dim_; i++ ) {
1348 if ( ((*ex)[i] >= x_up_[i]-epsn) ) {
1349 (*ev)[i] = 0.0;
1350 }
1351 }
1352 }
1353
1354 void pruneActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &x, Real eps = Real(0)) {
1355 ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1356 ROL::Ptr<std::vector<Real> > ev; cast_vector(ev,v);
1357 Real epsn = std::min(scale_*eps,min_diff_);
1358 for ( int i = 0; i < dim_; i++ ) {
1359 if ( ((*ex)[i] <= x_lo_[i]+epsn) ||
1360 ((*ex)[i] >= x_up_[i]-epsn) ) {
1361 (*ev)[i] = 0.0;
1362 }
1363 }
1364 }
1365
1366 void pruneLowerActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &g, const ROL::Vector<Real> &x, Real xeps = Real(0), Real geps = Real(0)) {
1367 ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1368 ROL::Ptr<const std::vector<Real> > eg; cast_const_vector(eg,g);
1369 ROL::Ptr<std::vector<Real> > ev; cast_vector(ev,v);
1370 Real epsn = std::min(scale_*xeps,min_diff_);
1371 for ( int i = 0; i < dim_; i++ ) {
1372 if ( ((*ex)[i] <= x_lo_[i]+epsn && (*eg)[i] > geps) ) {
1373 (*ev)[i] = 0.0;
1374 }
1375 }
1376 }
1377
1378 void pruneUpperActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &g, const ROL::Vector<Real> &x, Real xeps = Real(0), Real geps = Real(0)) {
1379 ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1380 ROL::Ptr<const std::vector<Real> > eg; cast_const_vector(eg,g);
1381 ROL::Ptr<std::vector<Real> > ev; cast_vector(ev,v);
1382 Real epsn = std::min(scale_*xeps,min_diff_);
1383 for ( int i = 0; i < dim_; i++ ) {
1384 if ( ((*ex)[i] >= x_up_[i]-epsn && (*eg)[i] < -geps) ) {
1385 (*ev)[i] = 0.0;
1386 }
1387 }
1388 }
1389
1390 void pruneActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &g, const ROL::Vector<Real> &x, Real xeps = Real(0), Real geps = Real(0)) {
1391 ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1392 ROL::Ptr<const std::vector<Real> > eg; cast_const_vector(eg,g);
1393 ROL::Ptr<std::vector<Real> > ev; cast_vector(ev,v);
1394 Real epsn = std::min(scale_*xeps,min_diff_);
1395 for ( int i = 0; i < dim_; i++ ) {
1396 if ( ((*ex)[i] <= x_lo_[i]+epsn && (*eg)[i] > geps) ||
1397 ((*ex)[i] >= x_up_[i]-epsn && (*eg)[i] < -geps) ) {
1398 (*ev)[i] = 0.0;
1399 }
1400 }
1401 }
1402
1403 const ROL::Ptr<const ROL::Vector<Real> > getLowerBound(void) const {
1404 return l_;
1405 }
1406
1407 const ROL::Ptr<const ROL::Vector<Real> > getUpperBound(void) const {
1408 return u_;
1409 }
1410};
1411
1412template<class Real>
1413class H1BoundConstraint : public ROL::BoundConstraint<Real> {
1414private:
1415 int dim_;
1416 std::vector<Real> x_lo_;
1417 std::vector<Real> x_up_;
1418 Real min_diff_;
1419 Real scale_;
1420 ROL::Ptr<BurgersFEM<Real> > fem_;
1421 ROL::Ptr<ROL::Vector<Real> > l_;
1422 ROL::Ptr<ROL::Vector<Real> > u_;
1423
1424 void cast_vector(ROL::Ptr<std::vector<Real> > &xvec,
1425 ROL::Vector<Real> &x) const {
1426 try {
1427 xvec = ROL::constPtrCast<std::vector<Real> >(
1428 (dynamic_cast<H1VectorPrimal<Real>&>(x)).getVector());
1429 }
1430 catch (std::exception &e) {
1431 xvec = ROL::constPtrCast<std::vector<Real> >(
1432 (dynamic_cast<H1VectorDual<Real>&>(x)).getVector());
1433 }
1434 }
1435
1436 void cast_const_vector(ROL::Ptr<const std::vector<Real> > &xvec,
1437 const ROL::Vector<Real> &x) const {
1438 try {
1439 xvec = (dynamic_cast<H1VectorPrimal<Real>&>(
1440 const_cast<ROL::Vector<Real> &>(x))).getVector();
1441 }
1442 catch (std::exception &e) {
1443 xvec = (dynamic_cast<H1VectorDual<Real>&>(
1444 const_cast<ROL::Vector<Real> &>(x))).getVector();
1445 }
1446 }
1447
1448 void axpy(std::vector<Real> &out, const Real a,
1449 const std::vector<Real> &x, const std::vector<Real> &y) const{
1450 out.resize(dim_,0.0);
1451 for (unsigned i = 0; i < dim_; i++) {
1452 out[i] = a*x[i] + y[i];
1453 }
1454 }
1455
1456 void projection(std::vector<Real> &x) {
1457 for ( int i = 0; i < dim_; i++ ) {
1458 x[i] = std::max(x_lo_[i],std::min(x_up_[i],x[i]));
1459 }
1460 }
1461
1462public:
1463 H1BoundConstraint(std::vector<Real> &l, std::vector<Real> &u,
1464 const ROL::Ptr<BurgersFEM<Real> > &fem, Real scale = 1.0)
1465 : x_lo_(l), x_up_(u), scale_(scale), fem_(fem) {
1466 dim_ = x_lo_.size();
1467 for ( int i = 0; i < dim_; i++ ) {
1468 if ( i == 0 ) {
1469 min_diff_ = x_up_[i] - x_lo_[i];
1470 }
1471 else {
1472 min_diff_ = ( (min_diff_ < (x_up_[i] - x_lo_[i])) ? min_diff_ : (x_up_[i] - x_lo_[i]) );
1473 }
1474 }
1475 min_diff_ *= 0.5;
1476 l_ = ROL::makePtr<H1VectorPrimal<Real>>(
1477 ROL::makePtr<std::vector<Real>>(l), fem);
1478 u_ = ROL::makePtr<H1VectorPrimal<Real>>(
1479 ROL::makePtr<std::vector<Real>>(u), fem);
1480 }
1481
1483 ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1484 bool val = true;
1485 int cnt = 1;
1486 for ( int i = 0; i < dim_; i++ ) {
1487 if ( (*ex)[i] >= x_lo_[i] && (*ex)[i] <= x_up_[i] ) { cnt *= 1; }
1488 else { cnt *= 0; }
1489 }
1490 if ( cnt == 0 ) { val = false; }
1491 return val;
1492 }
1493
1495 ROL::Ptr<std::vector<Real> > ex; cast_vector(ex,x);
1496 projection(*ex);
1497 }
1498
1499 void pruneLowerActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &x, Real eps = Real(0)) {
1500 ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1501 ROL::Ptr<std::vector<Real> > ev; cast_vector(ev,v);
1502 Real epsn = std::min(scale_*eps,min_diff_);
1503 for ( int i = 0; i < dim_; i++ ) {
1504 if ( ((*ex)[i] <= x_lo_[i]+epsn) ) {
1505 (*ev)[i] = 0.0;
1506 }
1507 }
1508 }
1509
1510 void pruneUpperActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &x, Real eps = Real(0)) {
1511 ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1512 ROL::Ptr<std::vector<Real> > ev; cast_vector(ev,v);
1513 Real epsn = std::min(scale_*eps,min_diff_);
1514 for ( int i = 0; i < dim_; i++ ) {
1515 if ( ((*ex)[i] >= x_up_[i]-epsn) ) {
1516 (*ev)[i] = 0.0;
1517 }
1518 }
1519 }
1520
1521 void pruneActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &x, Real eps = Real(0)) {
1522 ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1523 ROL::Ptr<std::vector<Real> > ev; cast_vector(ev,v);
1524 Real epsn = std::min(scale_*eps,min_diff_);
1525 for ( int i = 0; i < dim_; i++ ) {
1526 if ( ((*ex)[i] <= x_lo_[i]+epsn) ||
1527 ((*ex)[i] >= x_up_[i]-epsn) ) {
1528 (*ev)[i] = 0.0;
1529 }
1530 }
1531 }
1532
1533 void pruneLowerActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &g, const ROL::Vector<Real> &x, Real xeps = Real(0), Real geps = Real(0)) {
1534 ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1535 ROL::Ptr<const std::vector<Real> > eg; cast_const_vector(eg,g);
1536 ROL::Ptr<std::vector<Real> > ev; cast_vector(ev,v);
1537 Real epsn = std::min(scale_*xeps,min_diff_);
1538 for ( int i = 0; i < dim_; i++ ) {
1539 if ( ((*ex)[i] <= x_lo_[i]+epsn && (*eg)[i] > geps) ) {
1540 (*ev)[i] = 0.0;
1541 }
1542 }
1543 }
1544
1545 void pruneUpperActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &g, const ROL::Vector<Real> &x, Real xeps = Real(0), Real geps = Real(0)) {
1546 ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1547 ROL::Ptr<const std::vector<Real> > eg; cast_const_vector(eg,g);
1548 ROL::Ptr<std::vector<Real> > ev; cast_vector(ev,v);
1549 Real epsn = std::min(scale_*xeps,min_diff_);
1550 for ( int i = 0; i < dim_; i++ ) {
1551 if ( ((*ex)[i] >= x_up_[i]-epsn && (*eg)[i] < -geps) ) {
1552 (*ev)[i] = 0.0;
1553 }
1554 }
1555 }
1556
1557 void pruneActive(ROL::Vector<Real> &v, const ROL::Vector<Real> &g, const ROL::Vector<Real> &x, Real xeps = Real(0), Real geps = Real(0)) {
1558 ROL::Ptr<const std::vector<Real> > ex; cast_const_vector(ex,x);
1559 ROL::Ptr<const std::vector<Real> > eg; cast_const_vector(eg,g);
1560 ROL::Ptr<std::vector<Real> > ev; cast_vector(ev,v);
1561 Real epsn = std::min(scale_*xeps,min_diff_);
1562 for ( int i = 0; i < dim_; i++ ) {
1563 if ( ((*ex)[i] <= x_lo_[i]+epsn && (*eg)[i] > geps) ||
1564 ((*ex)[i] >= x_up_[i]-epsn && (*eg)[i] < -geps) ) {
1565 (*ev)[i] = 0.0;
1566 }
1567 }
1568 }
1569
1570 const ROL::Ptr<const ROL::Vector<Real> > getLowerBound(void) const {
1571 return l_;
1572 }
1573
1574 const ROL::Ptr<const ROL::Vector<Real> > getUpperBound(void) const {
1575 return u_;
1576 }
1577};
1578
1579template<class Real, class Ordinal>
1580class L2VectorBatchManager : public ROL::TeuchosBatchManager<Real,Ordinal> {
1581private:
1582 void cast_vector(ROL::Ptr<std::vector<Real> > &xvec,
1583 ROL::Vector<Real> &x) const {
1584 try {
1585 xvec = ROL::constPtrCast<std::vector<Real> >(
1586 (dynamic_cast<L2VectorPrimal<Real>&>(x)).getVector());
1587 }
1588 catch (std::exception &e) {
1589 xvec = ROL::constPtrCast<std::vector<Real> >(
1590 (dynamic_cast<L2VectorDual<Real>&>(x)).getVector());
1591 }
1592 }
1593
1594public:
1595 L2VectorBatchManager(const ROL::Ptr<const Teuchos::Comm<int> > &comm)
1596 : ROL::TeuchosBatchManager<Real,Ordinal>(comm) {}
1598 ROL::Ptr<std::vector<Real> > input_ptr;
1599 cast_vector(input_ptr,input);
1600 int dim_i = input_ptr->size();
1601 ROL::Ptr<std::vector<Real> > output_ptr;
1602 cast_vector(output_ptr,output);
1603 int dim_o = output_ptr->size();
1604 if ( dim_i != dim_o ) {
1605 std::cout << "L2VectorBatchManager: DIMENSION MISMATCH ON RANK "
1606 << ROL::TeuchosBatchManager<Real,Ordinal>::batchID() << "\n";
1607 }
1608 else {
1609 ROL::TeuchosBatchManager<Real,Ordinal>::sumAll(&(*input_ptr)[0],&(*output_ptr)[0],dim_i);
1610 }
1611 }
1612};
1613
1614template<class Real, class Ordinal>
1615class H1VectorBatchManager : public ROL::TeuchosBatchManager<Real,Ordinal> {
1616private:
1617 void cast_vector(ROL::Ptr<std::vector<Real> > &xvec,
1618 ROL::Vector<Real> &x) const {
1619 try {
1620 xvec = ROL::constPtrCast<std::vector<Real> >(
1621 (dynamic_cast<H1VectorPrimal<Real>&>(x)).getVector());
1622 }
1623 catch (std::exception &e) {
1624 xvec = ROL::constPtrCast<std::vector<Real> >(
1625 (dynamic_cast<H1VectorDual<Real>&>(x)).getVector());
1626 }
1627 }
1628
1629public:
1630 H1VectorBatchManager(const ROL::Ptr<const Teuchos::Comm<int> > &comm)
1631 : ROL::TeuchosBatchManager<Real,Ordinal>(comm) {}
1633 ROL::Ptr<std::vector<Real> > input_ptr;
1634 cast_vector(input_ptr,input);
1635 int dim_i = input_ptr->size();
1636 ROL::Ptr<std::vector<Real> > output_ptr;
1637 cast_vector(output_ptr,output);
1638 int dim_o = output_ptr->size();
1639 if ( dim_i != dim_o ) {
1640 std::cout << "H1VectorBatchManager: DIMENSION MISMATCH ON RANK "
1641 << ROL::TeuchosBatchManager<Real,Ordinal>::batchID() << "\n";
1642 }
1643 else {
1644 ROL::TeuchosBatchManager<Real,Ordinal>::sumAll(&(*input_ptr)[0],&(*output_ptr)[0],dim_i);
1645 }
1646 }
1647};
1648
1649template<class Real>
1650Real random(const ROL::Ptr<const Teuchos::Comm<int> > &comm) {
1651 Real val = 0.0;
1652 if ( Teuchos::rank<int>(*comm)==0 ) {
1653 val = (Real)rand()/(Real)RAND_MAX;
1654 }
1655 Teuchos::broadcast<int,Real>(*comm,0,1,&val);
1656 return val;
1657}
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
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)
Real get_viscosity(void) const
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)
BurgersFEM(int nx=128, Real nl=1.0, Real cH1=1.0, Real cL2=1.0)
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 set_problem_data(const Real nu, const Real f, const Real u0, const Real u1)
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 ...
std::vector< Real >::size_type uint
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.
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)
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
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 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 .
Real value(const ROL::Vector< Real > &z, Real &tol)
Compute value.
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.
Objective_BurgersControl(const ROL::Ptr< BurgersFEM< Real > > &fem, Real x=0.0)
std::vector< int > indices_
std::vector< Real >::size_type uint
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)
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.
const std::vector< Real > getParameter(void) const
Provides the interface to evaluate simulation-based objective functions.
Defines the linear algebra or vector space interface.
virtual void zero()
Set to zero vector.
Real random(const ROL::Ptr< const Teuchos::Comm< int > > &comm)
constexpr auto dim