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_02.hpp 45 \brief Provides definitions of equality constraint and objective for 46 example_02. 47 */ 48 49 #include "ROL_StdVector.hpp" 50 #include "ROL_Vector_SimOpt.hpp" 51 #include "ROL_Constraint_SimOpt.hpp" 52 #include "ROL_Objective_SimOpt.hpp" 53 #include "ROL_Reduced_Objective_SimOpt.hpp" 54 55 #include "ROL_Stream.hpp" 56 #include "Teuchos_GlobalMPISession.hpp" 57 #include "Teuchos_XMLParameterListHelpers.hpp" 58 #include "Teuchos_LAPACK.hpp" 59 60 #include <iostream> 61 #include <algorithm> 62 #include <ctime> 63 64 template<class Real> 65 class Constraint_BurgersControl : public ROL::Constraint_SimOpt<Real> { 66 private: 67 unsigned nx_; 68 unsigned nt_; 69 70 Real dx_; 71 Real T_; 72 Real dt_; 73 74 Real nu_; 75 Real u0_; 76 Real u1_; 77 Real f_; 78 std::vector<Real> u_init_; 79 80 private: compute_norm(const std::vector<Real> & r)81 Real compute_norm(const std::vector<Real> &r) { 82 return std::sqrt(dot(r,r)); 83 } 84 dot(const std::vector<Real> & x,const std::vector<Real> & y)85 Real dot(const std::vector<Real> &x, const std::vector<Real> &y) { 86 Real ip = 0.0; 87 Real c = ((x.size()==nx_) ? 4.0 : 2.0); 88 for (unsigned i = 0; i < x.size(); i++) { 89 if ( i == 0 ) { 90 ip += dx_/6.0*(c*x[i] + x[i+1])*y[i]; 91 } 92 else if ( i == x.size()-1 ) { 93 ip += dx_/6.0*(x[i-1] + c*x[i])*y[i]; 94 } 95 else { 96 ip += dx_/6.0*(x[i-1] + 4.0*x[i] + x[i+1])*y[i]; 97 } 98 } 99 return ip; 100 } 101 102 using ROL::Constraint_SimOpt<Real>::update; 103 update(std::vector<Real> & u,const std::vector<Real> & s,const Real alpha=1.0)104 void update(std::vector<Real> &u, const std::vector<Real> &s, const Real alpha=1.0) { 105 for (unsigned i = 0; i < u.size(); i++) { 106 u[i] += alpha*s[i]; 107 } 108 } 109 scale(std::vector<Real> & u,const Real alpha=0.0)110 void scale(std::vector<Real> &u, const Real alpha=0.0) { 111 for (unsigned i = 0; i < u.size(); i++) { 112 u[i] *= alpha; 113 } 114 } 115 compute_residual(std::vector<Real> & r,const std::vector<Real> & uold,const std::vector<Real> & zold,const std::vector<Real> & unew,const std::vector<Real> & znew)116 void compute_residual(std::vector<Real> &r, const std::vector<Real> &uold, const std::vector<Real> &zold, 117 const std::vector<Real> &unew, const std::vector<Real> &znew) { 118 r.clear(); 119 r.resize(nx_,0.0); 120 for (unsigned n = 0; n < nx_; n++) { 121 // Contribution from mass & stiffness term at time step t and t-1 122 r[n] += (4.0*dx_/6.0 + 0.5*dt_*2.0*nu_/dx_)*unew[n]; 123 r[n] += (-4.0*dx_/6.0 + 0.5*dt_*2.0*nu_/dx_)*uold[n]; 124 if ( n > 0 ) { 125 r[n] += (dx_/6.0 - 0.5*dt_*nu_/dx_)*unew[n-1]; 126 r[n] -= (dx_/6.0 + 0.5*dt_*nu_/dx_)*uold[n-1]; 127 } 128 if ( n < nx_-1 ) { 129 r[n] += (dx_/6.0 - 0.5*dt_*nu_/dx_)*unew[n+1]; 130 r[n] -= (dx_/6.0 + 0.5*dt_*nu_/dx_)*uold[n+1]; 131 } 132 // Contribution from nonlinear term 133 if ( n > 0 ) { 134 r[n] -= 0.5*dt_*unew[n-1]*(unew[n-1]+unew[n])/6.0; 135 r[n] -= 0.5*dt_*uold[n-1]*(uold[n-1]+uold[n])/6.0; 136 } 137 if ( n < nx_-1 ){ 138 r[n] += 0.5*dt_*unew[n+1]*(unew[n]+unew[n+1])/6.0; 139 r[n] += 0.5*dt_*uold[n+1]*(uold[n]+uold[n+1])/6.0; 140 } 141 // Contribution from control 142 r[n] -= 0.5*dt_*dx_/6.0*(znew[n]+4.0*znew[n+1]+znew[n+2]); 143 r[n] -= 0.5*dt_*dx_/6.0*(zold[n]+4.0*zold[n+1]+zold[n+2]); 144 // Contribution from right-hand side 145 r[n] -= dt_*dx_*f_; 146 } 147 // Contribution from Dirichlet boundary terms 148 r[0] -= dt_*(0.5*u0_*(unew[ 0] + uold[ 0])/6.0 + u0_*u0_/6.0 + nu_*u0_/dx_); 149 r[nx_-1] += dt_*(0.5*u1_*(unew[nx_-1] + uold[nx_-1])/6.0 + u1_*u1_/6.0 - nu_*u1_/dx_); 150 } 151 compute_pde_jacobian(std::vector<Real> & dl,std::vector<Real> & d,std::vector<Real> & du,const std::vector<Real> & u)152 void compute_pde_jacobian(std::vector<Real> &dl, std::vector<Real> &d, std::vector<Real> &du, 153 const std::vector<Real> &u) { 154 // Get Diagonal and Off-Diagonal Entries of linear PDE Jacobian 155 d.clear(); 156 d.resize(nx_,4.0*dx_/6.0 + 0.5*dt_*nu_*2.0/dx_); 157 dl.clear(); 158 dl.resize(nx_-1,dx_/6.0 - 0.5*dt_*nu_/dx_); 159 du.clear(); 160 du.resize(nx_-1,dx_/6.0 - 0.5*dt_*nu_/dx_); 161 // Contribution from nonlinearity 162 for (unsigned n = 0; n < nx_; n++) { 163 if ( n < nx_-1 ) { 164 dl[n] += 0.5*dt_*(-2.0*u[n]-u[n+1])/6.0; 165 d[n] += 0.5*dt_*u[n+1]/6.0; 166 } 167 if ( n > 0 ) { 168 d[n] -= 0.5*dt_*u[n-1]/6.0; 169 du[n-1] += 0.5*dt_*(u[n-1]+2.0*u[n])/6.0; 170 } 171 } 172 // Contribution from Dirichlet boundary conditions 173 d[0] -= 0.5*dt_*u0_/6.0; 174 d[nx_-1] += 0.5*dt_*u1_/6.0; 175 } 176 apply_pde_jacobian_new(std::vector<Real> & jv,const std::vector<Real> & v,const std::vector<Real> & u,bool adjoint=false)177 void apply_pde_jacobian_new(std::vector<Real> &jv, const std::vector<Real> &v, const std::vector<Real> &u, 178 bool adjoint = false) { 179 jv.clear(); 180 jv.resize(nx_,0.0); 181 // Fill Jacobian times a vector 182 for (unsigned n = 0; n < nx_; n++) { 183 jv[n] = (4.0*dx_/6.0 + 0.5*dt_*nu_/dx_*2.0)*v[n]; // Mass & Stiffness 184 if ( n > 0 ) { 185 jv[n] += dx_/6.0*v[n-1]-0.5*dt_*nu_/dx_*v[n-1]; // Mass & Stiffness 186 if ( adjoint ) { 187 jv[n] -= 0.5*dt_*(u[n-1]/6.0*v[n]-(u[n-1]+2.0*u[n])/6.0*v[n-1]); // Nonlinearity 188 } 189 else { 190 jv[n] -= 0.5*dt_*(u[n-1]/6.0*v[n]+(u[n]+2.0*u[n-1])/6.0*v[n-1]); // Nonlinearity 191 } 192 } 193 if ( n < nx_-1 ) { 194 jv[n] += dx_/6.0*v[n+1]-0.5*dt_*nu_/dx_*v[n+1]; // Mass & Stiffness 195 if ( adjoint ) { 196 jv[n] += 0.5*dt_*(u[n+1]/6.0*v[n]-(u[n+1]+2.0*u[n])/6.0*v[n+1]); // Nonlinearity 197 } 198 else { 199 jv[n] += 0.5*dt_*(u[n+1]/6.0*v[n]+(u[n]+2.0*u[n+1])/6.0*v[n+1]); // Nonlinearity 200 } 201 } 202 } 203 jv[0] -= 0.5*dt_*u0_/6.0*v[0]; // Nonlinearity 204 jv[nx_-1] += 0.5*dt_*u1_/6.0*v[nx_-1]; // Nonlinearity 205 } 206 apply_pde_jacobian_old(std::vector<Real> & jv,const std::vector<Real> & v,const std::vector<Real> & u,bool adjoint=false)207 void apply_pde_jacobian_old(std::vector<Real> &jv, const std::vector<Real> &v, const std::vector<Real> &u, 208 bool adjoint = false) { 209 jv.clear(); 210 jv.resize(nx_,0.0); 211 // Fill Jacobian times a vector 212 for (unsigned n = 0; n < nx_; n++) { 213 jv[n] = (-4.0*dx_/6.0 + 0.5*dt_*nu_/dx_*2.0)*v[n]; // Mass & Stiffness 214 if ( n > 0 ) { 215 jv[n] += -dx_/6.0*v[n-1]-0.5*dt_*nu_/dx_*v[n-1]; // Mass & Stiffness 216 if ( adjoint ) { 217 jv[n] -= 0.5*dt_*(u[n-1]/6.0*v[n]-(u[n-1]+2.0*u[n])/6.0*v[n-1]); // Nonlinearity 218 } 219 else { 220 jv[n] -= 0.5*dt_*(u[n-1]/6.0*v[n]+(u[n]+2.0*u[n-1])/6.0*v[n-1]); // Nonlinearity 221 } 222 } 223 if ( n < nx_-1 ) { 224 jv[n] += -dx_/6.0*v[n+1]-0.5*dt_*nu_/dx_*v[n+1]; // Mass & Stiffness 225 if ( adjoint ) { 226 jv[n] += 0.5*dt_*(u[n+1]/6.0*v[n]-(u[n+1]+2.0*u[n])/6.0*v[n+1]); // Nonlinearity 227 } 228 else { 229 jv[n] += 0.5*dt_*(u[n+1]/6.0*v[n]+(u[n]+2.0*u[n+1])/6.0*v[n+1]); // Nonlinearity 230 } 231 } 232 } 233 jv[0] -= 0.5*dt_*u0_/6.0*v[0]; // Nonlinearity 234 jv[nx_-1] += 0.5*dt_*u1_/6.0*v[nx_-1]; // Nonlinearity 235 } 236 apply_pde_jacobian(std::vector<Real> & jv,const std::vector<Real> & vold,const std::vector<Real> & uold,const std::vector<Real> & vnew,const std::vector<Real> unew,bool adjoint=false)237 void apply_pde_jacobian(std::vector<Real> &jv, const std::vector<Real> &vold, const std::vector<Real> &uold, 238 const std::vector<Real> &vnew, const std::vector<Real> unew, bool adjoint = false) { 239 jv.clear(); 240 jv.resize(nx_,0.0); 241 // Fill Jacobian times a vector 242 for (unsigned n = 0; n < nx_; n++) { 243 jv[n] += (4.0*dx_/6.0+0.5*dt_*nu_/dx_*2.0)*vnew[n]; // Mass & Stiffness 244 jv[n] -= (4.0*dx_/6.0-0.5*dt_*nu_/dx_*2.0)*vold[n]; // Mass & Stiffness 245 if ( n > 0 ) { 246 jv[n] += dx_/6.0*vnew[n-1]-0.5*dt_*nu_/dx_*vnew[n-1]; // Mass & Stiffness 247 jv[n] -= dx_/6.0*vold[n-1]+0.5*dt_*nu_/dx_*vold[n-1]; // Mass & Stiffness 248 if ( adjoint ) { 249 jv[n] -= 0.5*dt_*(unew[n-1]/6.0*vnew[n]-(unew[n-1]+2.0*unew[n])/6.0*vnew[n-1]); // Nonlinearity 250 jv[n] -= 0.5*dt_*(uold[n-1]/6.0*vold[n]-(uold[n-1]+2.0*uold[n])/6.0*vold[n-1]); // Nonlinearity 251 } 252 else { 253 jv[n] -= 0.5*dt_*(unew[n-1]/6.0*vnew[n]+(unew[n]+2.0*unew[n-1])/6.0*vnew[n-1]); // Nonlinearity 254 jv[n] -= 0.5*dt_*(uold[n-1]/6.0*vold[n]+(uold[n]+2.0*uold[n-1])/6.0*vold[n-1]); // Nonlinearity 255 } 256 } 257 if ( n < nx_-1 ) { 258 jv[n] += dx_/6.0*vnew[n+1]-0.5*dt_*nu_/dx_*vnew[n+1]; // Mass & Stiffness 259 jv[n] -= dx_/6.0*vold[n+1]+0.5*dt_*nu_/dx_*vold[n+1]; // Mass & Stiffness 260 if ( adjoint ) { 261 jv[n] += 0.5*dt_*(unew[n+1]/6.0*vnew[n]-(unew[n+1]+2.0*unew[n])/6.0*vnew[n+1]); // Nonlinearity 262 jv[n] += 0.5*dt_*(uold[n+1]/6.0*vold[n]-(uold[n+1]+2.0*uold[n])/6.0*vold[n+1]); // Nonlinearity 263 } 264 else { 265 jv[n] += 0.5*dt_*(unew[n+1]/6.0*vnew[n]+(unew[n]+2.0*unew[n+1])/6.0*vnew[n+1]); // Nonlinearity 266 jv[n] += 0.5*dt_*(uold[n+1]/6.0*vold[n]+(uold[n]+2.0*uold[n+1])/6.0*vold[n+1]); // Nonlinearity 267 } 268 } 269 } 270 jv[0] -= 0.5*dt_*u0_/6.0*vnew[0]; // Nonlinearity 271 jv[0] -= 0.5*dt_*u0_/6.0*vold[0]; // Nonlinearity 272 jv[nx_-1] += 0.5*dt_*u1_/6.0*vnew[nx_-1]; // Nonlinearity 273 jv[nx_-1] += 0.5*dt_*u1_/6.0*vold[nx_-1]; // Nonlinearity 274 } 275 apply_pde_hessian(std::vector<Real> & hv,const std::vector<Real> & wold,const std::vector<Real> & vold,const std::vector<Real> & wnew,const std::vector<Real> & vnew)276 void apply_pde_hessian(std::vector<Real> &hv, const std::vector<Real> &wold, const std::vector<Real> &vold, 277 const std::vector<Real> &wnew, const std::vector<Real> &vnew) { 278 hv.clear(); 279 hv.resize(nx_,0.0); 280 for (unsigned n = 0; n < nx_; n++) { 281 if ( n > 0 ) { 282 hv[n] += 0.5*dt_*((wnew[n-1]*(vnew[n-1]+2.0*vnew[n]) - wnew[n]*vnew[n-1])/6.0); 283 hv[n] += 0.5*dt_*((wold[n-1]*(vold[n-1]+2.0*vold[n]) - wold[n]*vold[n-1])/6.0); 284 } 285 if ( n < nx_-1 ){ 286 hv[n] += 0.5*dt_*((wnew[n]*vnew[n+1] - wnew[n+1]*(2.0*vnew[n]+vnew[n+1]))/6.0); 287 hv[n] += 0.5*dt_*((wold[n]*vold[n+1] - wold[n+1]*(2.0*vold[n]+vold[n+1]))/6.0); 288 } 289 } 290 } 291 apply_control_jacobian(std::vector<Real> & jv,const std::vector<Real> & v,bool adjoint=false)292 void apply_control_jacobian(std::vector<Real> &jv, const std::vector<Real> &v, bool adjoint = false) { 293 jv.clear(); 294 unsigned dim = ((adjoint == true) ? nx_+2 : nx_); 295 jv.resize(dim,0.0); 296 for (unsigned n = 0; n < dim; n++) { 297 if ( adjoint ) { 298 if ( n == 0 ) { 299 jv[n] = -0.5*dt_*dx_/6.0*v[n]; 300 } 301 else if ( n == 1 ) { 302 jv[n] = -0.5*dt_*dx_/6.0*(4.0*v[n-1]+v[n]); 303 } 304 else if ( n == nx_ ) { 305 jv[n] = -0.5*dt_*dx_/6.0*(4.0*v[n-1]+v[n-2]); 306 } 307 else if ( n == nx_+1 ) { 308 jv[n] = -0.5*dt_*dx_/6.0*v[n-2]; 309 } 310 else { 311 jv[n] = -0.5*dt_*dx_/6.0*(v[n-2]+4.0*v[n-1]+v[n]); 312 } 313 } 314 else { 315 jv[n] -= 0.5*dt_*dx_/6.0*(v[n]+4.0*v[n+1]+v[n+2]); 316 } 317 } 318 } 319 run_newton(std::vector<Real> & u,const std::vector<Real> & znew,const std::vector<Real> & uold,const std::vector<Real> & zold)320 void run_newton(std::vector<Real> &u, const std::vector<Real> &znew, 321 const std::vector<Real> &uold, const std::vector<Real> &zold) { 322 u.clear(); 323 u.resize(nx_,0.0); 324 // Compute residual and residual norm 325 std::vector<Real> r(nx_,0.0); 326 compute_residual(r,uold,zold,u,znew); 327 Real rnorm = compute_norm(r); 328 // Define tolerances 329 Real rtol = 1.e2*ROL::ROL_EPSILON<Real>(); 330 unsigned maxit = 500; 331 // Initialize Jacobian storage 332 std::vector<Real> d(nx_,0.0); 333 std::vector<Real> dl(nx_-1,0.0); 334 std::vector<Real> du(nx_-1,0.0); 335 // Iterate Newton's method 336 Real alpha = 1.0, tmp = 0.0; 337 std::vector<Real> s(nx_,0.0); 338 std::vector<Real> utmp(nx_,0.0); 339 for (unsigned i = 0; i < maxit; i++) { 340 //std::cout << i << " " << rnorm << "\n"; 341 // Get Jacobian 342 compute_pde_jacobian(dl,d,du,u); 343 // Solve Newton system 344 linear_solve(s,dl,d,du,r); 345 // Perform line search 346 tmp = rnorm; 347 alpha = 1.0; 348 utmp.assign(u.begin(),u.end()); 349 update(utmp,s,-alpha); 350 compute_residual(r,uold,zold,utmp,znew); 351 rnorm = compute_norm(r); 352 while ( rnorm > (1.0-1.e-4*alpha)*tmp && alpha > std::sqrt(ROL::ROL_EPSILON<Real>()) ) { 353 alpha *= 0.5; 354 utmp.assign(u.begin(),u.end()); 355 update(utmp,s,-alpha); 356 compute_residual(r,uold,zold,utmp,znew); 357 rnorm = compute_norm(r); 358 } 359 // Update iterate 360 u.assign(utmp.begin(),utmp.end()); 361 if ( rnorm < rtol ) { 362 break; 363 } 364 } 365 } 366 linear_solve(std::vector<Real> & u,const std::vector<Real> & dl,const std::vector<Real> & d,const std::vector<Real> & du,const std::vector<Real> & r,const bool transpose=false)367 void linear_solve(std::vector<Real> &u, 368 const std::vector<Real> &dl, const std::vector<Real> &d, const std::vector<Real> &du, 369 const std::vector<Real> &r, const bool transpose = false) { 370 bool useLAPACK = false; 371 if ( useLAPACK ) { // DIRECT SOLVE: USE LAPACK 372 u.assign(r.begin(),r.end()); 373 // Store matrix diagonal & off-diagonals. 374 std::vector<Real> Dl(dl); 375 std::vector<Real> Du(du); 376 std::vector<Real> D(d); 377 // Perform LDL factorization 378 Teuchos::LAPACK<int,Real> lp; 379 std::vector<Real> Du2(nx_-2,0.0); 380 std::vector<int> ipiv(nx_,0); 381 int info; 382 int ldb = nx_; 383 int nhrs = 1; 384 lp.GTTRF(nx_,&Dl[0],&D[0],&Du[0],&Du2[0],&ipiv[0],&info); 385 char trans = ((transpose == true) ? 'T' : 'N'); 386 lp.GTTRS(trans,nx_,nhrs,&Dl[0],&D[0],&Du[0],&Du2[0],&ipiv[0],&u[0],ldb,&info); 387 } 388 else { // ITERATIVE SOLVE: USE GAUSS-SEIDEL 389 u.clear(); 390 u.resize(nx_,0.0); 391 unsigned maxit = 100; 392 Real rtol = std::min(1.e-12,1.e-4*std::sqrt(dot(r,r))); 393 //Real rtol = 1e-6; 394 Real resid = 0.0; 395 Real rnorm = 10.0*rtol; 396 for (unsigned i = 0; i < maxit; i++) { 397 for (unsigned n = 0; n < nx_; n++) { 398 u[n] = r[n]/d[n]; 399 if ( n < nx_-1 ) { 400 u[n] -= ((transpose == false) ? du[n] : dl[n])*u[n+1]/d[n]; 401 } 402 if ( n > 0 ) { 403 u[n] -= ((transpose == false) ? dl[n-1] : du[n-1])*u[n-1]/d[n]; 404 } 405 } 406 // Compute Residual 407 rnorm = 0.0; 408 for (unsigned n = 0; n < nx_; n++) { 409 resid = r[n] - d[n]*u[n]; 410 if ( n < nx_-1 ) { 411 resid -= ((transpose == false) ? du[n] : dl[n])*u[n+1]; 412 } 413 if ( n > 0 ) { 414 resid -= ((transpose == false) ? dl[n-1] : du[n-1])*u[n-1]; 415 } 416 rnorm += resid*resid; 417 } 418 rnorm = std::sqrt(rnorm); 419 if ( rnorm < rtol ) { 420 //std::cout << "\ni = " << i+1 << " rnorm = " << rnorm << "\n"; 421 break; 422 } 423 } 424 } 425 } 426 427 public: 428 Constraint_BurgersControl(int nx=128,int nt=100,Real T=1,Real nu=1.e-2,Real u0=0.0,Real u1=0.0,Real f=0.0)429 Constraint_BurgersControl(int nx = 128, int nt = 100, Real T = 1, 430 Real nu = 1.e-2, Real u0 = 0.0, Real u1 = 0.0, Real f = 0.0) 431 : nx_((unsigned)nx), nt_((unsigned)nt), T_(T), nu_(nu), u0_(u0), u1_(u1), f_(f) { 432 dx_ = 1.0/((Real)nx+1.0); 433 dt_ = T/((Real)nt); 434 u_init_.clear(); 435 u_init_.resize(nx_,0.0); 436 Real x = 0.0; 437 for (unsigned n = 0; n < nx_; n++) { 438 x = (Real)(n+1)*dx_; 439 u_init_[n] = ((x <= 0.5) ? 1.0 : 0.0); 440 } 441 } 442 value(ROL::Vector<Real> & c,const ROL::Vector<Real> & u,const ROL::Vector<Real> & z,Real & tol)443 void value(ROL::Vector<Real> &c, const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) { 444 ROL::Ptr<std::vector<Real> > cp = 445 dynamic_cast<ROL::StdVector<Real>&>(c).getVector(); 446 ROL::Ptr<const std::vector<Real> > up = 447 dynamic_cast<const ROL::StdVector<Real>&>(u).getVector(); 448 ROL::Ptr<const std::vector<Real> > zp = 449 dynamic_cast<const ROL::StdVector<Real>&>(z).getVector(); 450 // Initialize storage 451 std::vector<Real> C(nx_,0.0); 452 std::vector<Real> uold(u_init_); 453 std::vector<Real> unew(u_init_); 454 std::vector<Real> zold(nx_+2,0.0); 455 std::vector<Real> znew(nx_+2,0.0); 456 // Copy initial control 457 for (unsigned n = 0; n < nx_+2; n++) { 458 zold[n] = (*zp)[n]; 459 } 460 // Evaluate residual 461 for (unsigned t = 0; t < nt_; t++) { 462 // Copy current state at time step t 463 for (unsigned n = 0; n < nx_; n++) { 464 unew[n] = (*up)[t*nx_+n]; 465 } 466 // Copy current control at time step t 467 for (unsigned n = 0; n < nx_+2; n++) { 468 znew[n] = (*zp)[(t+1)*(nx_+2)+n]; 469 } 470 // Compute residual at time step t 471 compute_residual(C,uold,zold,unew,znew); 472 // Store residual at time step t 473 for (unsigned n = 0; n < nx_; n++) { 474 (*cp)[t*nx_+n] = C[n]; 475 } 476 // Copy previous state and control at time step t+1 477 uold.assign(unew.begin(),unew.end()); 478 zold.assign(znew.begin(),znew.end()); 479 } 480 } 481 solve(ROL::Vector<Real> & c,ROL::Vector<Real> & u,const ROL::Vector<Real> & z,Real & tol)482 void solve(ROL::Vector<Real> &c, ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) { 483 ROL::Ptr<std::vector<Real> > up = 484 dynamic_cast<ROL::StdVector<Real>&>(u).getVector(); 485 up->assign(up->size(),z.norm()/up->size()); 486 ROL::Ptr<const std::vector<Real> > zp = 487 dynamic_cast<const ROL::StdVector<Real>&>(z).getVector(); 488 // Initialize storage 489 std::vector<Real> uold(u_init_); 490 std::vector<Real> unew(u_init_); 491 std::vector<Real> zold(nx_+2,0.0); 492 std::vector<Real> znew(nx_+2,0.0); 493 // Copy initial control 494 for (unsigned n = 0; n < nx_+2; n++) { 495 zold[n] = (*zp)[n]; 496 } 497 // Time step using the trapezoidal rule 498 for (unsigned t = 0; t < nt_; t++) { 499 // Copy current control at time step t 500 for (unsigned n = 0; n < nx_+2; n++) { 501 znew[n] = (*zp)[(t+1)*(nx_+2)+n]; 502 } 503 // Solve nonlinear system at time step t 504 run_newton(unew,znew,uold,zold); 505 // store state at time step t 506 for (unsigned n = 0; n < nx_; n++) { 507 (*up)[t*nx_+n] = unew[n]; 508 } 509 // Copy previous state and control at time step t+1 510 uold.assign(unew.begin(),unew.end()); 511 zold.assign(znew.begin(),znew.end()); 512 } 513 value(c,u,z,tol); 514 } 515 applyJacobian_1(ROL::Vector<Real> & jv,const ROL::Vector<Real> & v,const ROL::Vector<Real> & u,const ROL::Vector<Real> & z,Real & tol)516 void applyJacobian_1(ROL::Vector<Real> &jv, const ROL::Vector<Real> &v, const ROL::Vector<Real> &u, 517 const ROL::Vector<Real> &z, Real &tol) { 518 ROL::Ptr<std::vector<Real> > jvp = 519 dynamic_cast<ROL::StdVector<Real>&>(jv).getVector(); 520 ROL::Ptr<const std::vector<Real> > vp = 521 dynamic_cast<const ROL::StdVector<Real>&>(v).getVector(); 522 ROL::Ptr<const std::vector<Real> > up = 523 dynamic_cast<const ROL::StdVector<Real>&>(u).getVector(); 524 std::vector<Real> J(nx_,0.0); 525 std::vector<Real> vold(nx_,0.0); 526 std::vector<Real> vnew(nx_,0.0); 527 std::vector<Real> uold(nx_,0.0); 528 std::vector<Real> unew(nx_,0.0); 529 for (unsigned t = 0; t < nt_; t++) { 530 for (unsigned n = 0; n < nx_; n++) { 531 unew[n] = (*up)[t*nx_+n]; 532 vnew[n] = (*vp)[t*nx_+n]; 533 } 534 apply_pde_jacobian(J,vold,uold,vnew,unew); 535 for (unsigned n = 0; n < nx_; n++) { 536 (*jvp)[t*nx_+n] = J[n]; 537 } 538 vold.assign(vnew.begin(),vnew.end()); 539 uold.assign(unew.begin(),unew.end()); 540 } 541 } 542 applyJacobian_2(ROL::Vector<Real> & jv,const ROL::Vector<Real> & v,const ROL::Vector<Real> & u,const ROL::Vector<Real> & z,Real & tol)543 void applyJacobian_2(ROL::Vector<Real> &jv, const ROL::Vector<Real> &v, const ROL::Vector<Real> &u, 544 const ROL::Vector<Real> &z, Real &tol) { 545 ROL::Ptr<std::vector<Real> > jvp = 546 dynamic_cast<ROL::StdVector<Real>&>(jv).getVector(); 547 ROL::Ptr<const std::vector<Real> > vp = 548 dynamic_cast<const ROL::StdVector<Real>&>(v).getVector(); 549 ROL::Ptr<const std::vector<Real> > zp = 550 dynamic_cast<const ROL::StdVector<Real>&>(z).getVector(); 551 std::vector<Real> vnew(nx_+2,0.0); 552 std::vector<Real> vold(nx_+2,0.0); 553 std::vector<Real> jnew(nx_,0.0); 554 std::vector<Real> jold(nx_,0.0); 555 for (unsigned n = 0; n < nx_+2; n++) { 556 vold[n] = (*vp)[n]; 557 } 558 apply_control_jacobian(jold,vold); 559 for (unsigned t = 0; t < nt_; t++) { 560 for (unsigned n = 0; n < nx_+2; n++) { 561 vnew[n] = (*vp)[(t+1)*(nx_+2)+n]; 562 } 563 apply_control_jacobian(jnew,vnew); 564 for (unsigned n = 0; n < nx_; n++) { 565 // Contribution from control 566 (*jvp)[t*nx_+n] = jnew[n] + jold[n]; 567 } 568 jold.assign(jnew.begin(),jnew.end()); 569 } 570 } 571 applyInverseJacobian_1(ROL::Vector<Real> & ijv,const ROL::Vector<Real> & v,const ROL::Vector<Real> & u,const ROL::Vector<Real> & z,Real & tol)572 void applyInverseJacobian_1(ROL::Vector<Real> &ijv, const ROL::Vector<Real> &v, const ROL::Vector<Real> &u, 573 const ROL::Vector<Real> &z, Real &tol) { 574 ROL::Ptr<std::vector<Real> > ijvp = 575 dynamic_cast<ROL::StdVector<Real>&>(ijv).getVector(); 576 ROL::Ptr<const std::vector<Real> > vp = 577 dynamic_cast<const ROL::StdVector<Real>&>(v).getVector(); 578 ROL::Ptr<const std::vector<Real> > up = 579 dynamic_cast<const ROL::StdVector<Real>&>(u).getVector(); 580 std::vector<Real> J(nx_,0.0); 581 std::vector<Real> r(nx_,0.0); 582 std::vector<Real> s(nx_,0.0); 583 std::vector<Real> uold(nx_,0.0); 584 std::vector<Real> unew(nx_,0.0); 585 std::vector<Real> d(nx_,0.0); 586 std::vector<Real> dl(nx_-1,0.0); 587 std::vector<Real> du(nx_-1,0.0); 588 for (unsigned t = 0; t < nt_; t++) { 589 // Build rhs. 590 for (unsigned n = 0; n < nx_; n++) { 591 r[n] = (*vp)[t*nx_+n]; 592 unew[n] = (*up)[t*nx_+n]; 593 } 594 apply_pde_jacobian_old(J,s,uold); 595 update(r,J,-1.0); 596 // Build Jacobian. 597 compute_pde_jacobian(dl,d,du,unew); 598 // Solve linear system. 599 linear_solve(s,dl,d,du,r); 600 // Copy solution. 601 for (unsigned n = 0; n < nx_; n++) { 602 (*ijvp)[t*nx_+n] = s[n]; 603 } 604 uold.assign(unew.begin(),unew.end()); 605 } 606 } 607 applyAdjointJacobian_1(ROL::Vector<Real> & ajv,const ROL::Vector<Real> & v,const ROL::Vector<Real> & u,const ROL::Vector<Real> & z,Real & tol)608 void applyAdjointJacobian_1(ROL::Vector<Real> &ajv, const ROL::Vector<Real> &v, const ROL::Vector<Real> &u, 609 const ROL::Vector<Real> &z, Real &tol) { 610 ROL::Ptr<std::vector<Real> > jvp = 611 dynamic_cast<ROL::StdVector<Real>&>(ajv).getVector(); 612 ROL::Ptr<const std::vector<Real> > vp = 613 dynamic_cast<const ROL::StdVector<Real>&>(v).getVector(); 614 ROL::Ptr<const std::vector<Real> > up = 615 dynamic_cast<const ROL::StdVector<Real>&>(u).getVector(); 616 std::vector<Real> J(nx_,0.0); 617 std::vector<Real> vold(nx_,0.0); 618 std::vector<Real> vnew(nx_,0.0); 619 std::vector<Real> unew(nx_,0.0); 620 for (unsigned t = nt_; t > 0; t--) { 621 for (unsigned n = 0; n < nx_; n++) { 622 unew[n] = (*up)[(t-1)*nx_+n]; 623 vnew[n] = (*vp)[(t-1)*nx_+n]; 624 } 625 apply_pde_jacobian(J,vold,unew,vnew,unew,true); 626 for (unsigned n = 0; n < nx_; n++) { 627 (*jvp)[(t-1)*nx_+n] = J[n]; 628 } 629 vold.assign(vnew.begin(),vnew.end()); 630 } 631 } 632 applyAdjointJacobian_2(ROL::Vector<Real> & jv,const ROL::Vector<Real> & v,const ROL::Vector<Real> & u,const ROL::Vector<Real> & z,Real & tol)633 void applyAdjointJacobian_2(ROL::Vector<Real> &jv, const ROL::Vector<Real> &v, const ROL::Vector<Real> &u, 634 const ROL::Vector<Real> &z, Real &tol) { 635 ROL::Ptr<std::vector<Real> > jvp = 636 dynamic_cast<ROL::StdVector<Real>&>(jv).getVector(); 637 ROL::Ptr<const std::vector<Real> > vp = 638 dynamic_cast<const ROL::StdVector<Real>&>(v).getVector(); 639 ROL::Ptr<const std::vector<Real> > zp = 640 dynamic_cast<const ROL::StdVector<Real>&>(z).getVector(); 641 std::vector<Real> vnew(nx_,0.0); 642 std::vector<Real> vold(nx_,0.0); 643 std::vector<Real> jnew(nx_+2,0.0); 644 std::vector<Real> jold(nx_+2,0.0); 645 for (unsigned t = nt_+1; t > 0; t--) { 646 for (unsigned n = 0; n < nx_; n++) { 647 if ( t > 1 ) { 648 vnew[n] = (*vp)[(t-2)*nx_+n]; 649 } 650 else { 651 vnew[n] = 0.0; 652 } 653 } 654 apply_control_jacobian(jnew,vnew,true); 655 for (unsigned n = 0; n < nx_+2; n++) { 656 // Contribution from control 657 (*jvp)[(t-1)*(nx_+2)+n] = jnew[n] + jold[n]; 658 } 659 jold.assign(jnew.begin(),jnew.end()); 660 } 661 } 662 applyInverseAdjointJacobian_1(ROL::Vector<Real> & ijv,const ROL::Vector<Real> & v,const ROL::Vector<Real> & u,const ROL::Vector<Real> & z,Real & tol)663 void applyInverseAdjointJacobian_1(ROL::Vector<Real> &ijv, const ROL::Vector<Real> &v, 664 const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) { 665 ROL::Ptr<std::vector<Real> > ijvp = 666 dynamic_cast<ROL::StdVector<Real>&>(ijv).getVector(); 667 ROL::Ptr<const std::vector<Real> > vp = 668 dynamic_cast<const ROL::StdVector<Real>&>(v).getVector(); 669 ROL::Ptr<const std::vector<Real> > up = 670 dynamic_cast<const ROL::StdVector<Real>&>(u).getVector(); 671 std::vector<Real> J(nx_,0.0); 672 std::vector<Real> r(nx_,0.0); 673 std::vector<Real> s(nx_,0.0); 674 std::vector<Real> unew(nx_,0.0); 675 std::vector<Real> d(nx_,0.0); 676 std::vector<Real> dl(nx_-1,0.0); 677 std::vector<Real> du(nx_-1,0.0); 678 for (unsigned t = nt_; t > 0; t--) { 679 // Build rhs. 680 for (unsigned n = 0; n < nx_; n++) { 681 r[n] = (*vp)[(t-1)*nx_+n]; 682 unew[n] = (*up)[(t-1)*nx_+n]; 683 } 684 apply_pde_jacobian_old(J,s,unew,true); 685 update(r,J,-1.0); 686 // Build Jacobian. 687 compute_pde_jacobian(dl,d,du,unew); 688 // Solve linear system. 689 linear_solve(s,dl,d,du,r,true); 690 // Copy solution. 691 for (unsigned n = 0; n < nx_; n++) { 692 (*ijvp)[(t-1)*nx_+n] = s[n]; 693 } 694 } 695 } 696 applyAdjointHessian_11(ROL::Vector<Real> & hwv,const ROL::Vector<Real> & w,const ROL::Vector<Real> & v,const ROL::Vector<Real> & u,const ROL::Vector<Real> & z,Real & tol)697 void applyAdjointHessian_11(ROL::Vector<Real> &hwv, const ROL::Vector<Real> &w, const ROL::Vector<Real> &v, 698 const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) { 699 ROL::Ptr<std::vector<Real> > hwvp = 700 dynamic_cast<ROL::StdVector<Real>&>(hwv).getVector(); 701 ROL::Ptr<const std::vector<Real> > wp = 702 dynamic_cast<const ROL::StdVector<Real>&>(w).getVector(); 703 ROL::Ptr<const std::vector<Real> > vp = 704 dynamic_cast<const ROL::StdVector<Real>&>(v).getVector(); 705 std::vector<Real> snew(nx_,0.0); 706 std::vector<Real> wnew(nx_,0.0); 707 std::vector<Real> wold(nx_,0.0); 708 std::vector<Real> vnew(nx_,0.0); 709 for (unsigned t = nt_; t > 0; t--) { 710 for (unsigned n = 0; n < nx_; n++) { 711 vnew[n] = (*vp)[(t-1)*nx_+n]; 712 wnew[n] = (*wp)[(t-1)*nx_+n]; 713 } 714 apply_pde_hessian(snew,wold,vnew,wnew,vnew); 715 for (unsigned n = 0; n < nx_; n++) { 716 (*hwvp)[(t-1)*nx_+n] = snew[n]; 717 } 718 wold.assign(wnew.begin(),wnew.end()); 719 } 720 } 721 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)722 void applyAdjointHessian_12(ROL::Vector<Real> &ahwv, const ROL::Vector<Real> &w, const ROL::Vector<Real> &v, 723 const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) { 724 ahwv.zero(); 725 } 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)726 void applyAdjointHessian_21(ROL::Vector<Real> &ahwv, const ROL::Vector<Real> &w, const ROL::Vector<Real> &v, 727 const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) { 728 ahwv.zero(); 729 } 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)730 void applyAdjointHessian_22(ROL::Vector<Real> &ahwv, const ROL::Vector<Real> &w, const ROL::Vector<Real> &v, 731 const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) { 732 ahwv.zero(); 733 } 734 }; 735 736 template<class Real> 737 class Objective_BurgersControl : public ROL::Objective_SimOpt<Real> { 738 private: 739 Real alpha_; // Penalty Parameter 740 741 unsigned nx_; // Number of interior nodes 742 Real dx_; // Mesh spacing (i.e. 1/(nx+1)) 743 unsigned nt_; // Number of time steps 744 Real dt_; // Time step size 745 Real T_; // Final time 746 747 /***************************************************************/ 748 /********** BEGIN PRIVATE MEMBER FUNCTION DECLARATION **********/ 749 /***************************************************************/ evaluate_target(Real x)750 Real evaluate_target(Real x) { 751 Real val = 0.0; 752 int example = 1; 753 switch (example) { 754 case 1: val = ((x<=0.5) ? 1.0 : 0.0); break; 755 case 2: val = 1.0; break; 756 case 3: val = std::abs(std::sin(8.0*M_PI*x)); break; 757 case 4: val = std::exp(-0.5*(x-0.5)*(x-0.5)); break; 758 } 759 return val; 760 } 761 dot(const std::vector<Real> & x,const std::vector<Real> & y)762 Real dot(const std::vector<Real> &x, const std::vector<Real> &y) { 763 Real ip = 0.0; 764 Real c = ((x.size()==nx_) ? 4.0 : 2.0); 765 for (unsigned i=0; i<x.size(); i++) { 766 if ( i == 0 ) { 767 ip += dx_/6.0*(c*x[i] + x[i+1])*y[i]; 768 } 769 else if ( i == x.size()-1 ) { 770 ip += dx_/6.0*(x[i-1] + c*x[i])*y[i]; 771 } 772 else { 773 ip += dx_/6.0*(x[i-1] + 4.0*x[i] + x[i+1])*y[i]; 774 } 775 } 776 return ip; 777 } 778 apply_mass(std::vector<Real> & Mu,const std::vector<Real> & u)779 void apply_mass(std::vector<Real> &Mu, const std::vector<Real> &u ) { 780 Mu.resize(u.size(),0.0); 781 Real c = ((u.size()==nx_) ? 4.0 : 2.0); 782 for (unsigned i=0; i<u.size(); i++) { 783 if ( i == 0 ) { 784 Mu[i] = dx_/6.0*(c*u[i] + u[i+1]); 785 } 786 else if ( i == u.size()-1 ) { 787 Mu[i] = dx_/6.0*(u[i-1] + c*u[i]); 788 } 789 else { 790 Mu[i] = dx_/6.0*(u[i-1] + 4.0*u[i] + u[i+1]); 791 } 792 } 793 } 794 /*************************************************************/ 795 /********** END PRIVATE MEMBER FUNCTION DECLARATION **********/ 796 /*************************************************************/ 797 798 public: 799 Objective_BurgersControl(Real alpha=1.e-4,int nx=128,int nt=100,Real T=1.0)800 Objective_BurgersControl(Real alpha = 1.e-4, int nx = 128, int nt = 100, Real T = 1.0) 801 : alpha_(alpha), nx_((unsigned)nx), nt_((unsigned)nt), T_(T) { 802 dx_ = 1.0/((Real)nx+1.0); 803 dt_ = T/((Real)nt); 804 } 805 value(const ROL::Vector<Real> & u,const ROL::Vector<Real> & z,Real & tol)806 Real value( const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) { 807 ROL::Ptr<const std::vector<Real> > up = 808 dynamic_cast<const ROL::StdVector<Real>&>(u).getVector(); 809 ROL::Ptr<const std::vector<Real> > zp = 810 dynamic_cast<const ROL::StdVector<Real>&>(z).getVector(); 811 // COMPUTE RESIDUAL 812 std::vector<Real> U(nx_,0.0); 813 std::vector<Real> G(nx_,0.0); 814 std::vector<Real> Z(nx_+2,0.0); 815 for (unsigned n = 0; n < nx_+2; n++) { 816 Z[n] = (*zp)[n]; 817 } 818 Real ss = 0.5*dt_; 819 Real val = 0.5*ss*alpha_*dot(Z,Z); 820 for (unsigned t = 0; t < nt_; t++) { 821 ss = ((t < nt_-1) ? dt_ : 0.5*dt_); 822 for (unsigned n = 0; n < nx_; n++) { 823 U[n] = (*up)[t*nx_+n]-evaluate_target((Real)(n+1)*dx_); 824 G[n] = evaluate_target((Real)(n+1)*dx_); 825 } 826 val += 0.5*ss*dot(U,U); 827 val -= 0.5*ss*dot(G,G); // subtract constant term 828 for (unsigned n = 0; n < nx_+2; n++) { 829 Z[n] = (*zp)[(t+1)*(nx_+2)+n]; 830 } 831 val += 0.5*ss*alpha_*dot(Z,Z); 832 } 833 return val; 834 } 835 gradient_1(ROL::Vector<Real> & g,const ROL::Vector<Real> & u,const ROL::Vector<Real> & z,Real & tol)836 void gradient_1( ROL::Vector<Real> &g, const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) { 837 ROL::Ptr<std::vector<Real> > gp = 838 dynamic_cast<ROL::StdVector<Real>&>(g).getVector(); 839 ROL::Ptr<const std::vector<Real> > up = 840 dynamic_cast<const ROL::StdVector<Real>&>(u).getVector(); 841 // COMPUTE GRADIENT WRT U 842 std::vector<Real> U(nx_,0.0); 843 std::vector<Real> M(nx_,0.0); 844 Real ss = dt_; 845 for (unsigned t = 0; t < nt_; t++) { 846 ss = ((t < nt_-1) ? dt_ : 0.5*dt_); 847 for (unsigned n = 0; n < nx_; n++) { 848 U[n] = (*up)[t*nx_+n]-evaluate_target((Real)(n+1)*dx_); 849 } 850 apply_mass(M,U); 851 for (unsigned n = 0; n < nx_; n++) { 852 (*gp)[t*nx_+n] = ss*M[n]; 853 } 854 } 855 } 856 gradient_2(ROL::Vector<Real> & g,const ROL::Vector<Real> & u,const ROL::Vector<Real> & z,Real & tol)857 void gradient_2( ROL::Vector<Real> &g, const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) { 858 ROL::Ptr<std::vector<Real> > gp = 859 dynamic_cast<ROL::StdVector<Real>&>(g).getVector(); 860 ROL::Ptr<const std::vector<Real> > zp = 861 dynamic_cast<const ROL::StdVector<Real>&>(z).getVector(); 862 // COMPUTE GRADIENT WRT Z 863 std::vector<Real> Z(nx_+2,0.0); 864 std::vector<Real> M(nx_+2,0.0); 865 Real ss = dt_; 866 for (unsigned t = 0; t < nt_+1; t++) { 867 ss = ((t < nt_ && t > 0) ? dt_ : 0.5*dt_); 868 for (unsigned n = 0; n < nx_+2; n++) { 869 Z[n] = (*zp)[t*(nx_+2)+n]; 870 } 871 apply_mass(M,Z); 872 for (unsigned n = 0; n < nx_+2; n++) { 873 (*gp)[t*(nx_+2)+n] = ss*alpha_*M[n]; 874 } 875 } 876 } 877 hessVec_11(ROL::Vector<Real> & hv,const ROL::Vector<Real> & v,const ROL::Vector<Real> & u,const ROL::Vector<Real> & z,Real & tol)878 void hessVec_11( ROL::Vector<Real> &hv, const ROL::Vector<Real> &v, 879 const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) { 880 ROL::Ptr<std::vector<Real> > hvp = 881 dynamic_cast<ROL::StdVector<Real>&>(hv).getVector(); 882 ROL::Ptr<const std::vector<Real> > vp = 883 dynamic_cast<const ROL::StdVector<Real>&>(v).getVector(); 884 // COMPUTE GRADIENT WRT U 885 std::vector<Real> V(nx_,0.0); 886 std::vector<Real> M(nx_,0.0); 887 Real ss = 0.5*dt_; 888 for (unsigned t = 0; t < nt_; t++) { 889 ss = ((t < nt_-1) ? dt_ : 0.5*dt_); 890 for (unsigned n = 0; n < nx_; n++) { 891 V[n] = (*vp)[t*nx_+n]; 892 } 893 apply_mass(M,V); 894 for (unsigned n = 0; n < nx_; n++) { 895 (*hvp)[t*nx_+n] = ss*M[n]; 896 } 897 } 898 } 899 hessVec_12(ROL::Vector<Real> & hv,const ROL::Vector<Real> & v,const ROL::Vector<Real> & u,const ROL::Vector<Real> & z,Real & tol)900 void hessVec_12( ROL::Vector<Real> &hv, const ROL::Vector<Real> &v, 901 const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) { 902 hv.zero(); 903 } 904 hessVec_21(ROL::Vector<Real> & hv,const ROL::Vector<Real> & v,const ROL::Vector<Real> & u,const ROL::Vector<Real> & z,Real & tol)905 void hessVec_21( ROL::Vector<Real> &hv, const ROL::Vector<Real> &v, 906 const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) { 907 hv.zero(); 908 } 909 hessVec_22(ROL::Vector<Real> & hv,const ROL::Vector<Real> & v,const ROL::Vector<Real> & u,const ROL::Vector<Real> & z,Real & tol)910 void hessVec_22( ROL::Vector<Real> &hv, const ROL::Vector<Real> &v, 911 const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) { 912 ROL::Ptr<std::vector<Real> > hvp = ROL::constPtrCast<std::vector<Real> >( 913 (dynamic_cast<const ROL::StdVector<Real>&>(hv)).getVector()); 914 ROL::Ptr<const std::vector<Real> > vp = 915 dynamic_cast<const ROL::StdVector<Real>&>(v).getVector(); 916 // COMPUTE GRADIENT WRT Z 917 std::vector<Real> V(nx_+2,0.0); 918 std::vector<Real> M(nx_+2,0.0); 919 Real ss = 0.0; 920 for (unsigned t = 0; t < nt_+1; t++) { 921 ss = ((t < nt_ && t > 0) ? dt_ : 0.5*dt_); 922 for (unsigned n = 0; n < nx_+2; n++) { 923 V[n] = (*vp)[t*(nx_+2)+n]; 924 } 925 apply_mass(M,V); 926 for (unsigned n = 0; n < nx_+2; n++) { 927 (*hvp)[t*(nx_+2)+n] = ss*alpha_*M[n]; 928 } 929 } 930 } 931 }; 932