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