ROL
example_06.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 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::Vector<Real> {
532private:
533 ROL::Ptr<std::vector<Real> > vec_;
534 ROL::Ptr<BurgersFEM<Real> > fem_;
535
536 mutable ROL::Ptr<L2VectorDual<Real> > dual_vec_;
537
538public:
539 L2VectorPrimal(const ROL::Ptr<std::vector<Real> > & vec,
540 const ROL::Ptr<BurgersFEM<Real> > &fem)
541 : vec_(vec), fem_(fem), dual_vec_(ROL::nullPtr) {}
542
543 void set( const ROL::Vector<Real> &x ) {
544 const L2VectorPrimal &ex = dynamic_cast<const L2VectorPrimal&>(x);
545 const std::vector<Real>& xval = *ex.getVector();
546 std::copy(xval.begin(),xval.end(),vec_->begin());
547 }
548
549 void plus( const ROL::Vector<Real> &x ) {
550 const L2VectorPrimal &ex = dynamic_cast<const L2VectorPrimal&>(x);
551 const std::vector<Real>& xval = *ex.getVector();
552 unsigned dimension = vec_->size();
553 for (unsigned i=0; i<dimension; i++) {
554 (*vec_)[i] += xval[i];
555 }
556 }
557
558 void scale( const Real alpha ) {
559 unsigned dimension = vec_->size();
560 for (unsigned i=0; i<dimension; i++) {
561 (*vec_)[i] *= alpha;
562 }
563 }
564
565 Real dot( const ROL::Vector<Real> &x ) const {
566 const L2VectorPrimal & ex = dynamic_cast<const L2VectorPrimal&>(x);
567 const std::vector<Real>& xval = *ex.getVector();
568 return fem_->compute_L2_dot(xval,*vec_);
569 }
570
571 Real norm() const {
572 Real val = 0;
573 val = std::sqrt( dot(*this) );
574 return val;
575 }
576
577 ROL::Ptr<ROL::Vector<Real> > clone() const {
578 return ROL::makePtr<L2VectorPrimal>( ROL::makePtr<std::vector<Real>>(vec_->size(),0.0),fem_);
579 }
580
581 ROL::Ptr<const std::vector<Real> > getVector() const {
582 return vec_;
583 }
584
585 ROL::Ptr<std::vector<Real> > getVector() {
586 return vec_;
587 }
588
589 ROL::Ptr<ROL::Vector<Real> > basis( const int i ) const {
590 ROL::Ptr<L2VectorPrimal> e
591 = ROL::makePtr<L2VectorPrimal>( ROL::makePtr<std::vector<Real>>(vec_->size(),0.0),fem_);
592 (*e->getVector())[i] = 1.0;
593 return e;
594 }
595
596 int dimension() const {
597 return vec_->size();
598 }
599
600 const ROL::Vector<Real>& dual() const {
601 dual_vec_ = ROL::makePtr<L2VectorDual<Real>>(
602 ROL::makePtr<std::vector<Real>>(*vec_),fem_);
603
604 fem_->apply_mass(*(ROL::constPtrCast<std::vector<Real> >(dual_vec_->getVector())),*vec_);
605 return *dual_vec_;
606 }
607
608 Real apply(const ROL::Vector<Real> &x) const {
609 const L2VectorDual<Real> &ex = dynamic_cast<const L2VectorDual<Real>&>(x);
610 const std::vector<Real>& xval = *ex.getVector();
611 return std::inner_product(vec_->begin(), vec_->end(), xval.begin(), Real(0));
612 }
613
614};
615
616template<class Real>
617class L2VectorDual : public ROL::Vector<Real> {
618private:
619 ROL::Ptr<std::vector<Real> > vec_;
620 ROL::Ptr<BurgersFEM<Real> > fem_;
621
622 mutable ROL::Ptr<L2VectorPrimal<Real> > dual_vec_;
623
624public:
625 L2VectorDual(const ROL::Ptr<std::vector<Real> > & vec,
626 const ROL::Ptr<BurgersFEM<Real> > &fem)
627 : vec_(vec), fem_(fem), dual_vec_(ROL::nullPtr) {}
628
629 void set( const ROL::Vector<Real> &x ) {
630 const L2VectorDual &ex = dynamic_cast<const L2VectorDual&>(x);
631 const std::vector<Real>& xval = *ex.getVector();
632 std::copy(xval.begin(),xval.end(),vec_->begin());
633 }
634
635 void plus( const ROL::Vector<Real> &x ) {
636 const L2VectorDual &ex = dynamic_cast<const L2VectorDual&>(x);
637 const std::vector<Real>& xval = *ex.getVector();
638 unsigned dimension = vec_->size();
639 for (unsigned i=0; i<dimension; i++) {
640 (*vec_)[i] += xval[i];
641 }
642 }
643
644 void scale( const Real alpha ) {
645 unsigned dimension = vec_->size();
646 for (unsigned i=0; i<dimension; i++) {
647 (*vec_)[i] *= alpha;
648 }
649 }
650
651 Real dot( const ROL::Vector<Real> &x ) const {
652 const L2VectorDual & ex = dynamic_cast<const L2VectorDual&>(x);
653 const std::vector<Real>& xval = *ex.getVector();
654 unsigned dimension = vec_->size();
655 std::vector<Real> Mx(dimension,0.0);
656 fem_->apply_inverse_mass(Mx,xval);
657 Real val = 0.0;
658 for (unsigned i = 0; i < dimension; i++) {
659 val += Mx[i]*(*vec_)[i];
660 }
661 return val;
662 }
663
664 Real norm() const {
665 Real val = 0;
666 val = std::sqrt( dot(*this) );
667 return val;
668 }
669
670 ROL::Ptr<ROL::Vector<Real> > clone() const {
671 return ROL::makePtr<L2VectorDual>( ROL::makePtr<std::vector<Real>>(vec_->size(),0.0),fem_);
672 }
673
674 ROL::Ptr<const std::vector<Real> > getVector() const {
675 return vec_;
676 }
677
678 ROL::Ptr<std::vector<Real> > getVector() {
679 return vec_;
680 }
681
682 ROL::Ptr<ROL::Vector<Real> > basis( const int i ) const {
683 ROL::Ptr<L2VectorDual> e
684 = ROL::makePtr<L2VectorDual>( ROL::makePtr<std::vector<Real>>(vec_->size(),0.0),fem_);
685 (*e->getVector())[i] = 1.0;
686 return e;
687 }
688
689 int dimension() const {
690 return vec_->size();
691 }
692
693 const ROL::Vector<Real>& dual() const {
694 dual_vec_ = ROL::makePtr<L2VectorPrimal<Real>>(
695 ROL::makePtr<std::vector<Real>>(*vec_),fem_);
696
697 fem_->apply_inverse_mass(*(ROL::constPtrCast<std::vector<Real> >(dual_vec_->getVector())),*vec_);
698 return *dual_vec_;
699 }
700
701 Real apply(const ROL::Vector<Real> &x) const {
702 const L2VectorPrimal<Real> &ex = dynamic_cast<const L2VectorPrimal<Real>&>(x);
703 const std::vector<Real>& xval = *ex.getVector();
704 return std::inner_product(vec_->begin(), vec_->end(), xval.begin(), Real(0));
705 }
706
707};
708
709template<class Real>
710class H1VectorPrimal : public ROL::Vector<Real> {
711private:
712 ROL::Ptr<std::vector<Real> > vec_;
713 ROL::Ptr<BurgersFEM<Real> > fem_;
714
715 mutable ROL::Ptr<H1VectorDual<Real> > dual_vec_;
716
717public:
718 H1VectorPrimal(const ROL::Ptr<std::vector<Real> > & vec,
719 const ROL::Ptr<BurgersFEM<Real> > &fem)
720 : vec_(vec), fem_(fem), dual_vec_(ROL::nullPtr) {}
721
722 void set( const ROL::Vector<Real> &x ) {
723 const H1VectorPrimal &ex = dynamic_cast<const H1VectorPrimal&>(x);
724 const std::vector<Real>& xval = *ex.getVector();
725 std::copy(xval.begin(),xval.end(),vec_->begin());
726 }
727
728 void plus( const ROL::Vector<Real> &x ) {
729 const H1VectorPrimal &ex = dynamic_cast<const H1VectorPrimal&>(x);
730 const std::vector<Real>& xval = *ex.getVector();
731 unsigned dimension = vec_->size();
732 for (unsigned i=0; i<dimension; i++) {
733 (*vec_)[i] += xval[i];
734 }
735 }
736
737 void scale( const Real alpha ) {
738 unsigned dimension = vec_->size();
739 for (unsigned i=0; i<dimension; i++) {
740 (*vec_)[i] *= alpha;
741 }
742 }
743
744 Real dot( const ROL::Vector<Real> &x ) const {
745 const H1VectorPrimal & ex = dynamic_cast<const H1VectorPrimal&>(x);
746 const std::vector<Real>& xval = *ex.getVector();
747 return fem_->compute_H1_dot(xval,*vec_);
748 }
749
750 Real norm() const {
751 Real val = 0;
752 val = std::sqrt( dot(*this) );
753 return val;
754 }
755
756 ROL::Ptr<ROL::Vector<Real> > clone() const {
757 return ROL::makePtr<H1VectorPrimal>( ROL::makePtr<std::vector<Real>>(vec_->size(),0.0),fem_);
758 }
759
760 ROL::Ptr<const std::vector<Real> > getVector() const {
761 return vec_;
762 }
763
764 ROL::Ptr<std::vector<Real> > getVector() {
765 return vec_;
766 }
767
768 ROL::Ptr<ROL::Vector<Real> > basis( const int i ) const {
769 ROL::Ptr<H1VectorPrimal> e
770 = ROL::makePtr<H1VectorPrimal>( ROL::makePtr<std::vector<Real>>(vec_->size(),0.0),fem_);
771 (*e->getVector())[i] = 1.0;
772 return e;
773 }
774
775 int dimension() const {
776 return vec_->size();
777 }
778
779 const ROL::Vector<Real>& dual() const {
780 dual_vec_ = ROL::makePtr<H1VectorDual<Real>>(
781 ROL::makePtr<std::vector<Real>>(*vec_),fem_);
782
783 fem_->apply_H1(*(ROL::constPtrCast<std::vector<Real> >(dual_vec_->getVector())),*vec_);
784 return *dual_vec_;
785 }
786
787 Real apply(const ROL::Vector<Real> &x) const {
788 const H1VectorDual<Real> &ex = dynamic_cast<const H1VectorDual<Real>&>(x);
789 const std::vector<Real>& xval = *ex.getVector();
790 return std::inner_product(vec_->begin(), vec_->end(), xval.begin(), Real(0));
791 }
792
793};
794
795template<class Real>
796class H1VectorDual : public ROL::Vector<Real> {
797private:
798 ROL::Ptr<std::vector<Real> > vec_;
799 ROL::Ptr<BurgersFEM<Real> > fem_;
800
801 mutable ROL::Ptr<H1VectorPrimal<Real> > dual_vec_;
802
803public:
804 H1VectorDual(const ROL::Ptr<std::vector<Real> > & vec,
805 const ROL::Ptr<BurgersFEM<Real> > &fem)
806 : vec_(vec), fem_(fem), dual_vec_(ROL::nullPtr) {}
807
808 void set( const ROL::Vector<Real> &x ) {
809 const H1VectorDual &ex = dynamic_cast<const H1VectorDual&>(x);
810 const std::vector<Real>& xval = *ex.getVector();
811 std::copy(xval.begin(),xval.end(),vec_->begin());
812 }
813
814 void plus( const ROL::Vector<Real> &x ) {
815 const H1VectorDual &ex = dynamic_cast<const H1VectorDual&>(x);
816 const std::vector<Real>& xval = *ex.getVector();
817 unsigned dimension = vec_->size();
818 for (unsigned i=0; i<dimension; i++) {
819 (*vec_)[i] += xval[i];
820 }
821 }
822
823 void scale( const Real alpha ) {
824 unsigned dimension = vec_->size();
825 for (unsigned i=0; i<dimension; i++) {
826 (*vec_)[i] *= alpha;
827 }
828 }
829
830 Real dot( const ROL::Vector<Real> &x ) const {
831 const H1VectorDual & ex = dynamic_cast<const H1VectorDual&>(x);
832 const std::vector<Real>& xval = *ex.getVector();
833 unsigned dimension = vec_->size();
834 std::vector<Real> Mx(dimension,0.0);
835 fem_->apply_inverse_H1(Mx,xval);
836 Real val = 0.0;
837 for (unsigned i = 0; i < dimension; i++) {
838 val += Mx[i]*(*vec_)[i];
839 }
840 return val;
841 }
842
843 Real norm() const {
844 Real val = 0;
845 val = std::sqrt( dot(*this) );
846 return val;
847 }
848
849 ROL::Ptr<ROL::Vector<Real> > clone() const {
850 return ROL::makePtr<H1VectorDual>( ROL::makePtr<std::vector<Real>>(vec_->size(),0.0),fem_);
851 }
852
853 ROL::Ptr<const std::vector<Real> > getVector() const {
854 return vec_;
855 }
856
857 ROL::Ptr<std::vector<Real> > getVector() {
858 return vec_;
859 }
860
861 ROL::Ptr<ROL::Vector<Real> > basis( const int i ) const {
862 ROL::Ptr<H1VectorDual> e
863 = ROL::makePtr<H1VectorDual>( ROL::makePtr<std::vector<Real>>(vec_->size(),0.0),fem_);
864 (*e->getVector())[i] = 1.0;
865 return e;
866 }
867
868 int dimension() const {
869 return vec_->size();
870 }
871
872 const ROL::Vector<Real>& dual() const {
873 dual_vec_ = ROL::makePtr<H1VectorPrimal<Real>>(
874 ROL::makePtr<std::vector<Real>>(*vec_),fem_);
875
876 fem_->apply_inverse_H1(*(ROL::constPtrCast<std::vector<Real> >(dual_vec_->getVector())),*vec_);
877 return *dual_vec_;
878 }
879
880 Real apply(const ROL::Vector<Real> &x) const {
881 const H1VectorPrimal<Real> &ex = dynamic_cast<const H1VectorPrimal<Real>&>(x);
882 const std::vector<Real>& xval = *ex.getVector();
883 return std::inner_product(vec_->begin(), vec_->end(), xval.begin(), Real(0));
884 }
885
886};
887
888template<class Real>
890private:
891
894
897
900
901 ROL::Ptr<BurgersFEM<Real> > fem_;
902 bool useHessian_;
903
904public:
905 Constraint_BurgersControl(ROL::Ptr<BurgersFEM<Real> > &fem, bool useHessian = true)
906 : fem_(fem), useHessian_(useHessian) {}
907
909 const ROL::Vector<Real> &z, Real &tol) {
910 ROL::Ptr<std::vector<Real> > cp =
911 dynamic_cast<PrimalConstraintVector&>(c).getVector();
912 ROL::Ptr<const std::vector<Real> > up =
913 dynamic_cast<const PrimalStateVector&>(u).getVector();
914 ROL::Ptr<const std::vector<Real> > zp =
915 dynamic_cast<const PrimalControlVector&>(z).getVector();
916
917 const std::vector<Real> param
919 fem_->set_problem_data(param[0],param[1],param[2],param[3]);
920
921 fem_->compute_residual(*cp,*up,*zp);
922 }
923
925 const ROL::Vector<Real> &z, Real &tol) {
926 ROL::Ptr<std::vector<Real> > jvp =
927 dynamic_cast<PrimalConstraintVector&>(jv).getVector();
928 ROL::Ptr<const std::vector<Real> > vp =
929 dynamic_cast<const PrimalStateVector&>(v).getVector();
930 ROL::Ptr<const std::vector<Real> > up =
931 dynamic_cast<const PrimalStateVector&>(u).getVector();
932 ROL::Ptr<const std::vector<Real> > zp =
933 dynamic_cast<const PrimalControlVector&>(z).getVector();
934
935 const std::vector<Real> param
937 fem_->set_problem_data(param[0],param[1],param[2],param[3]);
938
939 fem_->apply_pde_jacobian(*jvp,*vp,*up,*zp);
940 }
941
943 const ROL::Vector<Real> &z, Real &tol) {
944 ROL::Ptr<std::vector<Real> > jvp =
945 dynamic_cast<PrimalConstraintVector&>(jv).getVector();
946 ROL::Ptr<const std::vector<Real> > vp =
947 dynamic_cast<const PrimalControlVector&>(v).getVector();
948 ROL::Ptr<const std::vector<Real> > up =
949 dynamic_cast<const PrimalStateVector&>(u).getVector();
950 ROL::Ptr<const std::vector<Real> > zp =
951 dynamic_cast<const PrimalControlVector&>(z).getVector();
952
953 const std::vector<Real> param
955 fem_->set_problem_data(param[0],param[1],param[2],param[3]);
956
957 fem_->apply_control_jacobian(*jvp,*vp,*up,*zp);
958 }
959
961 const ROL::Vector<Real> &z, Real &tol) {
962 ROL::Ptr<std::vector<Real> > ijvp =
963 dynamic_cast<PrimalStateVector&>(ijv).getVector();
964 ROL::Ptr<const std::vector<Real> > vp =
965 dynamic_cast<const PrimalConstraintVector&>(v).getVector();
966 ROL::Ptr<const std::vector<Real> > up =
967 dynamic_cast<const PrimalStateVector&>(u).getVector();
968 ROL::Ptr<const std::vector<Real> > zp =
969 dynamic_cast<const PrimalControlVector&>(z).getVector();
970
971 const std::vector<Real> param
973 fem_->set_problem_data(param[0],param[1],param[2],param[3]);
974
975 fem_->apply_inverse_pde_jacobian(*ijvp,*vp,*up,*zp);
976 }
977
979 const ROL::Vector<Real> &z, Real &tol) {
980 ROL::Ptr<std::vector<Real> > jvp =
981 dynamic_cast<DualStateVector&>(ajv).getVector();
982 ROL::Ptr<const std::vector<Real> > vp =
983 dynamic_cast<const DualConstraintVector&>(v).getVector();
984 ROL::Ptr<const std::vector<Real> > up =
985 dynamic_cast<const PrimalStateVector&>(u).getVector();
986 ROL::Ptr<const std::vector<Real> > zp =
987 dynamic_cast<const PrimalControlVector&>(z).getVector();
988
989 const std::vector<Real> param
991 fem_->set_problem_data(param[0],param[1],param[2],param[3]);
992
993 fem_->apply_adjoint_pde_jacobian(*jvp,*vp,*up,*zp);
994 }
995
997 const ROL::Vector<Real> &z, Real &tol) {
998 ROL::Ptr<std::vector<Real> > jvp =
999 dynamic_cast<DualControlVector&>(jv).getVector();
1000 ROL::Ptr<const std::vector<Real> > vp =
1001 dynamic_cast<const DualConstraintVector&>(v).getVector();
1002 ROL::Ptr<const std::vector<Real> > up =
1003 dynamic_cast<const PrimalStateVector&>(u).getVector();
1004 ROL::Ptr<const std::vector<Real> > zp =
1005 dynamic_cast<const PrimalControlVector&>(z).getVector();
1006
1007 const std::vector<Real> param
1009 fem_->set_problem_data(param[0],param[1],param[2],param[3]);
1010
1011 fem_->apply_adjoint_control_jacobian(*jvp,*vp,*up,*zp);
1012 }
1013
1015 const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
1016 ROL::Ptr<std::vector<Real> > iajvp =
1017 dynamic_cast<DualConstraintVector&>(iajv).getVector();
1018 ROL::Ptr<const std::vector<Real> > vp =
1019 dynamic_cast<const DualStateVector&>(v).getVector();
1020 ROL::Ptr<const std::vector<Real> > up =
1021 dynamic_cast<const PrimalStateVector&>(u).getVector();
1022 ROL::Ptr<const std::vector<Real> > zp =
1023 dynamic_cast<const PrimalControlVector&>(z).getVector();
1024
1025 const std::vector<Real> param
1027 fem_->set_problem_data(param[0],param[1],param[2],param[3]);
1028
1029 fem_->apply_inverse_adjoint_pde_jacobian(*iajvp,*vp,*up,*zp);
1030 }
1031
1033 const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
1034 if ( useHessian_ ) {
1035 ROL::Ptr<std::vector<Real> > ahwvp =
1036 dynamic_cast<DualStateVector&>(ahwv).getVector();
1037 ROL::Ptr<const std::vector<Real> > wp =
1038 dynamic_cast<const DualConstraintVector&>(w).getVector();
1039 ROL::Ptr<const std::vector<Real> > vp =
1040 dynamic_cast<const PrimalStateVector&>(v).getVector();
1041 ROL::Ptr<const std::vector<Real> > up =
1042 dynamic_cast<const PrimalStateVector&>(u).getVector();
1043 ROL::Ptr<const std::vector<Real> > zp =
1044 dynamic_cast<const PrimalControlVector&>(z).getVector();
1045
1046 const std::vector<Real> param
1048 fem_->set_problem_data(param[0],param[1],param[2],param[3]);
1049
1050 fem_->apply_adjoint_pde_hessian(*ahwvp,*wp,*vp,*up,*zp);
1051 }
1052 else {
1053 ahwv.zero();
1054 }
1055 }
1056
1058 const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
1059 if ( useHessian_ ) {
1060 ROL::Ptr<std::vector<Real> > ahwvp =
1061 dynamic_cast<DualControlVector&>(ahwv).getVector();
1062 ROL::Ptr<const std::vector<Real> > wp =
1063 dynamic_cast<const DualConstraintVector&>(w).getVector();
1064 ROL::Ptr<const std::vector<Real> > vp =
1065 dynamic_cast<const PrimalStateVector&>(v).getVector();
1066 ROL::Ptr<const std::vector<Real> > up =
1067 dynamic_cast<const PrimalStateVector&>(u).getVector();
1068 ROL::Ptr<const std::vector<Real> > zp =
1069 dynamic_cast<const PrimalControlVector&>(z).getVector();
1070
1071 const std::vector<Real> param
1073 fem_->set_problem_data(param[0],param[1],param[2],param[3]);
1074
1075 fem_->apply_adjoint_control_pde_hessian(*ahwvp,*wp,*vp,*up,*zp);
1076 }
1077 else {
1078 ahwv.zero();
1079 }
1080 }
1082 const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
1083 if ( useHessian_ ) {
1084 ROL::Ptr<std::vector<Real> > ahwvp =
1085 dynamic_cast<DualStateVector&>(ahwv).getVector();
1086 ROL::Ptr<const std::vector<Real> > wp =
1087 dynamic_cast<const DualConstraintVector&>(w).getVector();
1088 ROL::Ptr<const std::vector<Real> > vp =
1089 dynamic_cast<const PrimalControlVector&>(v).getVector();
1090 ROL::Ptr<const std::vector<Real> > up =
1091 dynamic_cast<const PrimalStateVector&>(u).getVector();
1092 ROL::Ptr<const std::vector<Real> > zp =
1093 dynamic_cast<const PrimalControlVector&>(z).getVector();
1094
1095 const std::vector<Real> param
1097 fem_->set_problem_data(param[0],param[1],param[2],param[3]);
1098
1099 fem_->apply_adjoint_pde_control_hessian(*ahwvp,*wp,*vp,*up,*zp);
1100 }
1101 else {
1102 ahwv.zero();
1103 }
1104 }
1106 const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
1107 if ( useHessian_ ) {
1108 ROL::Ptr<std::vector<Real> > ahwvp =
1109 dynamic_cast<DualControlVector&>(ahwv).getVector();
1110 ROL::Ptr<const std::vector<Real> > wp =
1111 dynamic_cast<const DualConstraintVector&>(w).getVector();
1112 ROL::Ptr<const std::vector<Real> > vp =
1113 dynamic_cast<const PrimalControlVector&>(v).getVector();
1114 ROL::Ptr<const std::vector<Real> > up =
1115 dynamic_cast<const PrimalStateVector&>(u).getVector();
1116 ROL::Ptr<const std::vector<Real> > zp =
1117 dynamic_cast<const PrimalControlVector&>(z).getVector();
1118
1119 const std::vector<Real> param
1121 fem_->set_problem_data(param[0],param[1],param[2],param[3]);
1122
1123 fem_->apply_adjoint_control_hessian(*ahwvp,*wp,*vp,*up,*zp);
1124 }
1125 else {
1126 ahwv.zero();
1127 }
1128 }
1129};
1130
1131template<class Real>
1133private:
1134
1137
1140
1141 Real alpha_; // Penalty Parameter
1142 ROL::Ptr<BurgersFEM<Real> > fem_;
1143 ROL::Ptr<ROL::Vector<Real> > ud_;
1144 ROL::Ptr<ROL::Vector<Real> > diff_;
1145
1146public:
1148 const ROL::Ptr<ROL::Vector<Real> > &ud,
1149 Real alpha = 1.e-4) : alpha_(alpha), fem_(fem), ud_(ud) {
1150 diff_ = ud_->clone();
1151 }
1152
1153 using ROL::Objective_SimOpt<Real>::value;
1154
1155 Real value( const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
1156 ROL::Ptr<const std::vector<Real> > up =
1157 dynamic_cast<const PrimalStateVector&>(u).getVector();
1158 ROL::Ptr<const std::vector<Real> > zp =
1159 dynamic_cast<const PrimalControlVector&>(z).getVector();
1160 ROL::Ptr<const std::vector<Real> > udp =
1161 dynamic_cast<const L2VectorPrimal<Real>&>(*ud_).getVector();
1162
1163 std::vector<Real> diff(udp->size(),0.0);
1164 for (unsigned i = 0; i < udp->size(); i++) {
1165 diff[i] = (*up)[i] - (*udp)[i];
1166 }
1167 return 0.5*(fem_->compute_L2_dot(diff,diff) + alpha_*fem_->compute_L2_dot(*zp,*zp));
1168 }
1169
1170 void gradient_1( ROL::Vector<Real> &g, const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
1171 ROL::Ptr<std::vector<Real> > gp =
1172 dynamic_cast<DualStateVector&>(g).getVector();
1173 ROL::Ptr<const std::vector<Real> > up =
1174 dynamic_cast<const PrimalStateVector&>(u).getVector();
1175 ROL::Ptr<const std::vector<Real> > udp =
1176 dynamic_cast<const L2VectorPrimal<Real>&>(*ud_).getVector();
1177
1178 std::vector<Real> diff(udp->size(),0.0);
1179 for (unsigned i = 0; i < udp->size(); i++) {
1180 diff[i] = (*up)[i] - (*udp)[i];
1181 }
1182 fem_->apply_mass(*gp,diff);
1183 }
1184
1185 void gradient_2( ROL::Vector<Real> &g, const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
1186 ROL::Ptr<std::vector<Real> > gp =
1187 dynamic_cast<DualControlVector&>(g).getVector();
1188 ROL::Ptr<const std::vector<Real> > zp =
1189 dynamic_cast<const PrimalControlVector&>(z).getVector();
1190
1191 fem_->apply_mass(*gp,*zp);
1192 for (unsigned i = 0; i < zp->size(); i++) {
1193 (*gp)[i] *= alpha_;
1194 }
1195 }
1196
1198 const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
1199 ROL::Ptr<std::vector<Real> > hvp =
1200 dynamic_cast<DualStateVector&>(hv).getVector();
1201 ROL::Ptr<const std::vector<Real> > vp =
1202 dynamic_cast<const PrimalStateVector&>(v).getVector();
1203
1204 fem_->apply_mass(*hvp,*vp);
1205 }
1206
1208 const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
1209 hv.zero();
1210 }
1211
1213 const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
1214 hv.zero();
1215 }
1216
1218 const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
1219 ROL::Ptr<std::vector<Real> > hvp =
1220 dynamic_cast<DualControlVector&>(hv).getVector();
1221 ROL::Ptr<const std::vector<Real> > vp =
1222 dynamic_cast<const PrimalControlVector&>(v).getVector();
1223
1224 fem_->apply_mass(*hvp,*vp);
1225 for (unsigned i = 0; i < vp->size(); i++) {
1226 (*hvp)[i] *= alpha_;
1227 }
1228 }
1229};
1230
1231template<class Real, class Ordinal>
1232class L2VectorBatchManager : public ROL::TeuchosBatchManager<Real,Ordinal> {
1233private:
1234 void cast_vector(ROL::Ptr<std::vector<Real> > &xvec,
1235 ROL::Vector<Real> &x) const {
1236 try {
1237 xvec = dynamic_cast<L2VectorPrimal<Real>&>(x).getVector();
1238 }
1239 catch (std::exception &e) {
1240 xvec = dynamic_cast<L2VectorDual<Real>&>(x).getVector();
1241 }
1242 }
1243
1244public:
1245 L2VectorBatchManager(const ROL::Ptr<const Teuchos::Comm<int> > &comm)
1246 : ROL::TeuchosBatchManager<Real,Ordinal>(comm) {}
1248 ROL::Ptr<std::vector<Real> > input_ptr;
1249 cast_vector(input_ptr,input);
1250 int dim_i = input_ptr->size();
1251 ROL::Ptr<std::vector<Real> > output_ptr;
1252 cast_vector(output_ptr,output);
1253 int dim_o = output_ptr->size();
1254 if ( dim_i != dim_o ) {
1255 std::cout << "L2VectorBatchManager: DIMENSION MISMATCH ON RANK "
1256 << ROL::TeuchosBatchManager<Real,Ordinal>::batchID() << "\n";
1257 }
1258 else {
1259 ROL::TeuchosBatchManager<Real,Ordinal>::sumAll(&(*input_ptr)[0],&(*output_ptr)[0],dim_i);
1260 }
1261 }
1262};
1263
1264template<class Real, class Ordinal>
1265class H1VectorBatchManager : public ROL::TeuchosBatchManager<Real,Ordinal> {
1266private:
1267 void cast_vector(ROL::Ptr<std::vector<Real> > &xvec,
1268 ROL::Vector<Real> &x) const {
1269 try {
1270 xvec = dynamic_cast<H1VectorPrimal<Real>&>(x).getVector();
1271 }
1272 catch (std::exception &e) {
1273 xvec = dynamic_cast<H1VectorDual<Real>&>(x).getVector();
1274 }
1275 }
1276
1277public:
1278 H1VectorBatchManager(const ROL::Ptr<const Teuchos::Comm<int> > &comm)
1279 : ROL::TeuchosBatchManager<Real,Ordinal>(comm) {}
1281 ROL::Ptr<std::vector<Real> > input_ptr;
1282 cast_vector(input_ptr,input);
1283 int dim_i = input_ptr->size();
1284 ROL::Ptr<std::vector<Real> > output_ptr;
1285 cast_vector(output_ptr,output);
1286 int dim_o = output_ptr->size();
1287 if ( dim_i != dim_o ) {
1288 std::cout << "H1VectorBatchManager: DIMENSION MISMATCH ON RANK "
1289 << ROL::TeuchosBatchManager<Real,Ordinal>::batchID() << "\n";
1290 }
1291 else {
1292 ROL::TeuchosBatchManager<Real,Ordinal>::sumAll(&(*input_ptr)[0],&(*output_ptr)[0],dim_i);
1293 }
1294 }
1295};
1296
1297template<class Real>
1298Real random(const ROL::Ptr<const Teuchos::Comm<int> > &comm) {
1299 Real val = 0.0;
1300 if ( Teuchos::rank<int>(*comm)==0 ) {
1301 val = (Real)rand()/(Real)RAND_MAX;
1302 }
1303 Teuchos::broadcast<int,Real>(*comm,0,1,&val);
1304 return val;
1305}
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)
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 ...
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 .
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 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 .
Objective_BurgersControl(const ROL::Ptr< BurgersFEM< Real > > &fem, const ROL::Ptr< ROL::Vector< Real > > &ud, Real alpha=1.e-4)
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.
ROL::Ptr< ROL::Vector< Real > > ud_
void hessVec_11(ROL::Vector< Real > &hv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Apply Hessian approximation to vector.
void hessVec_22(ROL::Vector< Real > &hv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
ROL::Ptr< ROL::Vector< Real > > diff_
void gradient_2(ROL::Vector< Real > &g, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Compute gradient with respect to second component.
L2VectorPrimal< Real > PrimalControlVector
ROL::Ptr< BurgersFEM< Real > > fem_
H1VectorDual< Real > DualStateVector
void hessVec_12(ROL::Vector< Real > &hv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Real value(const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Compute value.
L2VectorDual< Real > DualControlVector
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