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