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