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"
50 #include "ROL_Types.hpp"
51 
52 #include "ROL_StdVector.hpp"
53 #include "ROL_Vector_SimOpt.hpp"
55 #include "ROL_Objective_SimOpt.hpp"
57 
58 #include "Teuchos_oblackholestream.hpp"
59 #include "Teuchos_GlobalMPISession.hpp"
60 #include "Teuchos_XMLParameterListHelpers.hpp"
61 #include "Teuchos_LAPACK.hpp"
62 
63 #include <iostream>
64 #include <algorithm>
65 #include <ctime>
66 
67 template<class Real>
69 private:
70  unsigned nx_;
71  unsigned nt_;
72 
73  Real dx_;
74  Real T_;
75  Real dt_;
76 
77  Real nu_;
78  Real u0_;
79  Real u1_;
80  Real f_;
81  std::vector<Real> u_init_;
82 
83 private:
84  Real compute_norm(const std::vector<Real> &r) {
85  return std::sqrt(dot(r,r));
86  }
87 
88  Real dot(const std::vector<Real> &x, const std::vector<Real> &y) {
89  Real ip = 0.0;
90  Real c = ((x.size()==nx_) ? 4.0 : 2.0);
91  for (unsigned i = 0; i < x.size(); i++) {
92  if ( i == 0 ) {
93  ip += dx_/6.0*(c*x[i] + x[i+1])*y[i];
94  }
95  else if ( i == x.size()-1 ) {
96  ip += dx_/6.0*(x[i-1] + c*x[i])*y[i];
97  }
98  else {
99  ip += dx_/6.0*(x[i-1] + 4.0*x[i] + x[i+1])*y[i];
100  }
101  }
102  return ip;
103  }
104 
105  void update(std::vector<Real> &u, const std::vector<Real> &s, const Real alpha=1.0) {
106  for (unsigned i = 0; i < u.size(); i++) {
107  u[i] += alpha*s[i];
108  }
109  }
110 
111  void scale(std::vector<Real> &u, const Real alpha=0.0) {
112  for (unsigned i = 0; i < u.size(); i++) {
113  u[i] *= alpha;
114  }
115  }
116 
117  void compute_residual(std::vector<Real> &r, const std::vector<Real> &uold, const std::vector<Real> &zold,
118  const std::vector<Real> &unew, const std::vector<Real> &znew) {
119  r.clear();
120  r.resize(nx_,0.0);
121  for (unsigned n = 0; n < nx_; n++) {
122  // Contribution from mass & stiffness term at time step t and t-1
123  r[n] += (4.0*dx_/6.0 + 0.5*dt_*2.0*nu_/dx_)*unew[n];
124  r[n] += (-4.0*dx_/6.0 + 0.5*dt_*2.0*nu_/dx_)*uold[n];
125  if ( n > 0 ) {
126  r[n] += (dx_/6.0 - 0.5*dt_*nu_/dx_)*unew[n-1];
127  r[n] -= (dx_/6.0 + 0.5*dt_*nu_/dx_)*uold[n-1];
128  }
129  if ( n < nx_-1 ) {
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  // Contribution from nonlinear term
134  if ( n > 0 ) {
135  r[n] -= 0.5*dt_*unew[n-1]*(unew[n-1]+unew[n])/6.0;
136  r[n] -= 0.5*dt_*uold[n-1]*(uold[n-1]+uold[n])/6.0;
137  }
138  if ( n < nx_-1 ){
139  r[n] += 0.5*dt_*unew[n+1]*(unew[n]+unew[n+1])/6.0;
140  r[n] += 0.5*dt_*uold[n+1]*(uold[n]+uold[n+1])/6.0;
141  }
142  // Contribution from control
143  r[n] -= 0.5*dt_*dx_/6.0*(znew[n]+4.0*znew[n+1]+znew[n+2]);
144  r[n] -= 0.5*dt_*dx_/6.0*(zold[n]+4.0*zold[n+1]+zold[n+2]);
145  // Contribution from right-hand side
146  r[n] -= dt_*dx_*f_;
147  }
148  // Contribution from Dirichlet boundary terms
149  r[0] -= dt_*(0.5*u0_*(unew[ 0] + uold[ 0])/6.0 + u0_*u0_/6.0 + nu_*u0_/dx_);
150  r[nx_-1] += dt_*(0.5*u1_*(unew[nx_-1] + uold[nx_-1])/6.0 + u1_*u1_/6.0 - nu_*u1_/dx_);
151  }
152 
153  void compute_pde_jacobian(std::vector<Real> &dl, std::vector<Real> &d, std::vector<Real> &du,
154  const std::vector<Real> &u) {
155  // Get Diagonal and Off-Diagonal Entries of linear PDE Jacobian
156  d.clear();
157  d.resize(nx_,4.0*dx_/6.0 + 0.5*dt_*nu_*2.0/dx_);
158  dl.clear();
159  dl.resize(nx_-1,dx_/6.0 - 0.5*dt_*nu_/dx_);
160  du.clear();
161  du.resize(nx_-1,dx_/6.0 - 0.5*dt_*nu_/dx_);
162  // Contribution from nonlinearity
163  for (unsigned n = 0; n < nx_; n++) {
164  if ( n < nx_-1 ) {
165  dl[n] += 0.5*dt_*(-2.0*u[n]-u[n+1])/6.0;
166  d[n] += 0.5*dt_*u[n+1]/6.0;
167  }
168  if ( n > 0 ) {
169  d[n] -= 0.5*dt_*u[n-1]/6.0;
170  du[n-1] += 0.5*dt_*(u[n-1]+2.0*u[n])/6.0;
171  }
172  }
173  // Contribution from Dirichlet boundary conditions
174  d[0] -= 0.5*dt_*u0_/6.0;
175  d[nx_-1] += 0.5*dt_*u1_/6.0;
176  }
177 
178  void apply_pde_jacobian_new(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_old(std::vector<Real> &jv, const std::vector<Real> &v, const std::vector<Real> &u,
209  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)*v[n]; // Mass & Stiffness
215  if ( n > 0 ) {
216  jv[n] += -dx_/6.0*v[n-1]-0.5*dt_*nu_/dx_*v[n-1]; // Mass & Stiffness
217  if ( adjoint ) {
218  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
219  }
220  else {
221  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
222  }
223  }
224  if ( n < nx_-1 ) {
225  jv[n] += -dx_/6.0*v[n+1]-0.5*dt_*nu_/dx_*v[n+1]; // Mass & Stiffness
226  if ( adjoint ) {
227  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
228  }
229  else {
230  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
231  }
232  }
233  }
234  jv[0] -= 0.5*dt_*u0_/6.0*v[0]; // Nonlinearity
235  jv[nx_-1] += 0.5*dt_*u1_/6.0*v[nx_-1]; // Nonlinearity
236  }
237 
238  void apply_pde_jacobian(std::vector<Real> &jv, const std::vector<Real> &vold, const std::vector<Real> &uold,
239  const std::vector<Real> &vnew, const std::vector<Real> unew, bool adjoint = false) {
240  jv.clear();
241  jv.resize(nx_,0.0);
242  // Fill Jacobian times a vector
243  for (unsigned n = 0; n < nx_; n++) {
244  jv[n] += (4.0*dx_/6.0+0.5*dt_*nu_/dx_*2.0)*vnew[n]; // Mass & Stiffness
245  jv[n] -= (4.0*dx_/6.0-0.5*dt_*nu_/dx_*2.0)*vold[n]; // Mass & Stiffness
246  if ( n > 0 ) {
247  jv[n] += dx_/6.0*vnew[n-1]-0.5*dt_*nu_/dx_*vnew[n-1]; // Mass & Stiffness
248  jv[n] -= dx_/6.0*vold[n-1]+0.5*dt_*nu_/dx_*vold[n-1]; // Mass & Stiffness
249  if ( adjoint ) {
250  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
251  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
252  }
253  else {
254  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
255  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
256  }
257  }
258  if ( n < nx_-1 ) {
259  jv[n] += dx_/6.0*vnew[n+1]-0.5*dt_*nu_/dx_*vnew[n+1]; // Mass & Stiffness
260  jv[n] -= dx_/6.0*vold[n+1]+0.5*dt_*nu_/dx_*vold[n+1]; // Mass & Stiffness
261  if ( adjoint ) {
262  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
263  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
264  }
265  else {
266  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
267  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
268  }
269  }
270  }
271  jv[0] -= 0.5*dt_*u0_/6.0*vnew[0]; // Nonlinearity
272  jv[0] -= 0.5*dt_*u0_/6.0*vold[0]; // Nonlinearity
273  jv[nx_-1] += 0.5*dt_*u1_/6.0*vnew[nx_-1]; // Nonlinearity
274  jv[nx_-1] += 0.5*dt_*u1_/6.0*vold[nx_-1]; // Nonlinearity
275  }
276 
277  void apply_pde_hessian(std::vector<Real> &hv, const std::vector<Real> &wold, const std::vector<Real> &vold,
278  const std::vector<Real> &wnew, const std::vector<Real> &vnew) {
279  hv.clear();
280  hv.resize(nx_,0.0);
281  for (unsigned n = 0; n < nx_; n++) {
282  if ( n > 0 ) {
283  hv[n] += 0.5*dt_*((wnew[n-1]*(vnew[n-1]+2.0*vnew[n]) - wnew[n]*vnew[n-1])/6.0);
284  hv[n] += 0.5*dt_*((wold[n-1]*(vold[n-1]+2.0*vold[n]) - wold[n]*vold[n-1])/6.0);
285  }
286  if ( n < nx_-1 ){
287  hv[n] += 0.5*dt_*((wnew[n]*vnew[n+1] - wnew[n+1]*(2.0*vnew[n]+vnew[n+1]))/6.0);
288  hv[n] += 0.5*dt_*((wold[n]*vold[n+1] - wold[n+1]*(2.0*vold[n]+vold[n+1]))/6.0);
289  }
290  }
291  }
292 
293  void apply_control_jacobian(std::vector<Real> &jv, const std::vector<Real> &v, bool adjoint = false) {
294  jv.clear();
295  unsigned dim = ((adjoint == true) ? nx_+2 : nx_);
296  jv.resize(dim,0.0);
297  for (unsigned n = 0; n < dim; n++) {
298  if ( adjoint ) {
299  if ( n == 0 ) {
300  jv[n] = -0.5*dt_*dx_/6.0*v[n];
301  }
302  else if ( n == 1 ) {
303  jv[n] = -0.5*dt_*dx_/6.0*(4.0*v[n-1]+v[n]);
304  }
305  else if ( n == nx_ ) {
306  jv[n] = -0.5*dt_*dx_/6.0*(4.0*v[n-1]+v[n-2]);
307  }
308  else if ( n == nx_+1 ) {
309  jv[n] = -0.5*dt_*dx_/6.0*v[n-2];
310  }
311  else {
312  jv[n] = -0.5*dt_*dx_/6.0*(v[n-2]+4.0*v[n-1]+v[n]);
313  }
314  }
315  else {
316  jv[n] -= 0.5*dt_*dx_/6.0*(v[n]+4.0*v[n+1]+v[n+2]);
317  }
318  }
319  }
320 
321  void run_newton(std::vector<Real> &u, const std::vector<Real> &znew,
322  const std::vector<Real> &uold, const std::vector<Real> &zold) {
323  u.clear();
324  u.resize(nx_,0.0);
325  // Compute residual and residual norm
326  std::vector<Real> r(nx_,0.0);
327  compute_residual(r,uold,zold,u,znew);
328  Real rnorm = compute_norm(r);
329  // Define tolerances
330  Real rtol = 1.e2*ROL::ROL_EPSILON;
331  unsigned maxit = 500;
332  // Initialize Jacobian storage
333  std::vector<Real> d(nx_,0.0);
334  std::vector<Real> dl(nx_-1,0.0);
335  std::vector<Real> du(nx_-1,0.0);
336  // Iterate Newton's method
337  Real alpha = 1.0, tmp = 0.0;
338  std::vector<Real> s(nx_,0.0);
339  std::vector<Real> utmp(nx_,0.0);
340  for (unsigned i = 0; i < maxit; i++) {
341  //std::cout << i << " " << rnorm << "\n";
342  // Get Jacobian
343  compute_pde_jacobian(dl,d,du,u);
344  // Solve Newton system
345  linear_solve(s,dl,d,du,r);
346  // Perform line search
347  tmp = rnorm;
348  alpha = 1.0;
349  utmp.assign(u.begin(),u.end());
350  update(utmp,s,-alpha);
351  compute_residual(r,uold,zold,utmp,znew);
352  rnorm = compute_norm(r);
353  while ( rnorm > (1.0-1.e-4*alpha)*tmp && alpha > std::sqrt(ROL::ROL_EPSILON) ) {
354  alpha *= 0.5;
355  utmp.assign(u.begin(),u.end());
356  update(utmp,s,-alpha);
357  compute_residual(r,uold,zold,utmp,znew);
358  rnorm = compute_norm(r);
359  }
360  // Update iterate
361  u.assign(utmp.begin(),utmp.end());
362  if ( rnorm < rtol ) {
363  break;
364  }
365  }
366  }
367 
368  void linear_solve(std::vector<Real> &u,
369  const std::vector<Real> &dl, const std::vector<Real> &d, const std::vector<Real> &du,
370  const std::vector<Real> &r, const bool transpose = false) {
371  bool useLAPACK = false;
372  if ( useLAPACK ) { // DIRECT SOLVE: USE LAPACK
373  u.assign(r.begin(),r.end());
374  // Store matrix diagonal & off-diagonals.
375  std::vector<Real> Dl(dl);
376  std::vector<Real> Du(du);
377  std::vector<Real> D(d);
378  // Perform LDL factorization
379  Teuchos::LAPACK<int,Real> lp;
380  std::vector<Real> Du2(nx_-2,0.0);
381  std::vector<int> ipiv(nx_,0);
382  int info;
383  int ldb = nx_;
384  int nhrs = 1;
385  lp.GTTRF(nx_,&Dl[0],&D[0],&Du[0],&Du2[0],&ipiv[0],&info);
386  char trans = ((transpose == true) ? 'T' : 'N');
387  lp.GTTRS(trans,nx_,nhrs,&Dl[0],&D[0],&Du[0],&Du2[0],&ipiv[0],&u[0],ldb,&info);
388  }
389  else { // ITERATIVE SOLVE: USE GAUSS-SEIDEL
390  u.clear();
391  u.resize(nx_,0.0);
392  unsigned maxit = 100;
393  Real rtol = std::min(1.e-12,1.e-4*std::sqrt(dot(r,r)));
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 
429  EqualityConstraint_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 
443  void value(ROL::Vector<Real> &c, const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
444  Teuchos::RCP<std::vector<Real> > cp =
445  Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<ROL::StdVector<Real> >(c)).getVector());
446  Teuchos::RCP<const std::vector<Real> > up =
447  (Teuchos::dyn_cast<ROL::StdVector<Real> >(const_cast<ROL::Vector<Real> &>(u))).getVector();
448  Teuchos::RCP<const std::vector<Real> > zp =
449  (Teuchos::dyn_cast<ROL::StdVector<Real> >(const_cast<ROL::Vector<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 
482  void solve(ROL::Vector<Real> &c, ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
483  Teuchos::RCP<std::vector<Real> > up =
484  Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<ROL::StdVector<Real> >(u)).getVector());
485  up->assign(up->size(),z.norm()/up->size());
486  Teuchos::RCP<const std::vector<Real> > zp =
487  (Teuchos::dyn_cast<ROL::StdVector<Real> >(const_cast<ROL::Vector<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 
517  const ROL::Vector<Real> &z, Real &tol) {
518  Teuchos::RCP<std::vector<Real> > jvp =
519  Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<ROL::StdVector<Real> >(jv)).getVector());
520  Teuchos::RCP<const std::vector<Real> > vp =
521  (Teuchos::dyn_cast<ROL::StdVector<Real> >(const_cast<ROL::Vector<Real> &>(v))).getVector();
522  Teuchos::RCP<const std::vector<Real> > up =
523  (Teuchos::dyn_cast<ROL::StdVector<Real> >(const_cast<ROL::Vector<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 
544  const ROL::Vector<Real> &z, Real &tol) {
545  Teuchos::RCP<std::vector<Real> > jvp =
546  Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<ROL::StdVector<Real> >(jv)).getVector());
547  Teuchos::RCP<const std::vector<Real> > vp =
548  (Teuchos::dyn_cast<ROL::StdVector<Real> >(const_cast<ROL::Vector<Real> &>(v))).getVector();
549  Teuchos::RCP<const std::vector<Real> > zp =
550  (Teuchos::dyn_cast<ROL::StdVector<Real> >(const_cast<ROL::Vector<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 
573  const ROL::Vector<Real> &z, Real &tol) {
574  Teuchos::RCP<std::vector<Real> > ijvp =
575  Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<ROL::StdVector<Real> >(ijv)).getVector());
576  Teuchos::RCP<const std::vector<Real> > vp =
577  (Teuchos::dyn_cast<ROL::StdVector<Real> >(const_cast<ROL::Vector<Real> &>(v))).getVector();
578  Teuchos::RCP<const std::vector<Real> > up =
579  (Teuchos::dyn_cast<ROL::StdVector<Real> >(const_cast<ROL::Vector<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 
609  const ROL::Vector<Real> &z, Real &tol) {
610  Teuchos::RCP<std::vector<Real> > jvp =
611  Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<ROL::StdVector<Real> >(ajv)).getVector());
612  Teuchos::RCP<const std::vector<Real> > vp =
613  (Teuchos::dyn_cast<ROL::StdVector<Real> >(const_cast<ROL::Vector<Real> &>(v))).getVector();
614  Teuchos::RCP<const std::vector<Real> > up =
615  (Teuchos::dyn_cast<ROL::StdVector<Real> >(const_cast<ROL::Vector<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 
634  const ROL::Vector<Real> &z, Real &tol) {
635  Teuchos::RCP<std::vector<Real> > jvp =
636  Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<ROL::StdVector<Real> >(jv)).getVector());
637  Teuchos::RCP<const std::vector<Real> > vp =
638  (Teuchos::dyn_cast<ROL::StdVector<Real> >(const_cast<ROL::Vector<Real> &>(v))).getVector();
639  Teuchos::RCP<const std::vector<Real> > zp =
640  (Teuchos::dyn_cast<ROL::StdVector<Real> >(const_cast<ROL::Vector<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 
664  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
665  Teuchos::RCP<std::vector<Real> > ijvp =
666  Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<ROL::StdVector<Real> >(ijv)).getVector());
667  Teuchos::RCP<const std::vector<Real> > vp =
668  (Teuchos::dyn_cast<ROL::StdVector<Real> >(const_cast<ROL::Vector<Real> &>(v))).getVector();
669  Teuchos::RCP<const std::vector<Real> > up =
670  (Teuchos::dyn_cast<ROL::StdVector<Real> >(const_cast<ROL::Vector<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 
698  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
699  Teuchos::RCP<std::vector<Real> > hwvp =
700  Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<ROL::StdVector<Real> >(hwv)).getVector());
701  Teuchos::RCP<const std::vector<Real> > wp =
702  (Teuchos::dyn_cast<ROL::StdVector<Real> >(const_cast<ROL::Vector<Real> &>(w))).getVector();
703  Teuchos::RCP<const std::vector<Real> > vp =
704  (Teuchos::dyn_cast<ROL::StdVector<Real> >(const_cast<ROL::Vector<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 
723  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
724  ahwv.zero();
725  }
727  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
728  ahwv.zero();
729  }
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 /***************************************************************/
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 
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 
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 
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 
806  Real value( const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
807  Teuchos::RCP<const std::vector<Real> > up =
808  (Teuchos::dyn_cast<ROL::StdVector<Real> >(const_cast<ROL::Vector<Real> &>(u))).getVector();
809  Teuchos::RCP<const std::vector<Real> > zp =
810  (Teuchos::dyn_cast<ROL::StdVector<Real> >(const_cast<ROL::Vector<Real> &>(z))).getVector();
811  // COMPUTE RESIDUAL
812  std::vector<Real> U(nx_,0.0);
813  std::vector<Real> Z(nx_+2,0.0);
814  for (unsigned n = 0; n < nx_+2; n++) {
815  Z[n] = (*zp)[n];
816  }
817  Real ss = 0.5*dt_;
818  Real val = 0.5*ss*alpha_*dot(Z,Z);
819  for (unsigned t = 0; t < nt_; t++) {
820  ss = ((t < nt_-1) ? dt_ : 0.5*dt_);
821  for (unsigned n = 0; n < nx_; n++) {
822  U[n] = (*up)[t*nx_+n]-evaluate_target((Real)(n+1)*dx_);
823  }
824  val += 0.5*ss*dot(U,U);
825  for (unsigned n = 0; n < nx_+2; n++) {
826  Z[n] = (*zp)[(t+1)*(nx_+2)+n];
827  }
828  val += 0.5*ss*alpha_*dot(Z,Z);
829  }
830  return val;
831  }
832 
833  void gradient_1( ROL::Vector<Real> &g, const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
834  Teuchos::RCP<std::vector<Real> > gp = Teuchos::rcp_const_cast<std::vector<Real> >(
835  (Teuchos::dyn_cast<const ROL::StdVector<Real> >(g)).getVector());
836  Teuchos::RCP<const std::vector<Real> > up =
837  (Teuchos::dyn_cast<ROL::StdVector<Real> >(const_cast<ROL::Vector<Real> &>(u))).getVector();
838  // COMPUTE GRADIENT WRT U
839  std::vector<Real> U(nx_,0.0);
840  std::vector<Real> M(nx_,0.0);
841  Real ss = dt_;
842  for (unsigned t = 0; t < nt_; t++) {
843  ss = ((t < nt_-1) ? dt_ : 0.5*dt_);
844  for (unsigned n = 0; n < nx_; n++) {
845  U[n] = (*up)[t*nx_+n]-evaluate_target((Real)(n+1)*dx_);
846  }
847  apply_mass(M,U);
848  for (unsigned n = 0; n < nx_; n++) {
849  (*gp)[t*nx_+n] = ss*M[n];
850  }
851  }
852  }
853 
854  void gradient_2( ROL::Vector<Real> &g, const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
855  Teuchos::RCP<std::vector<Real> > gp = Teuchos::rcp_const_cast<std::vector<Real> >(
856  (Teuchos::dyn_cast<const ROL::StdVector<Real> >(g)).getVector());
857  Teuchos::RCP<const std::vector<Real> > zp =
858  (Teuchos::dyn_cast<ROL::StdVector<Real> >(const_cast<ROL::Vector<Real> &>(z))).getVector();
859  // COMPUTE GRADIENT WRT Z
860  std::vector<Real> Z(nx_+2,0.0);
861  std::vector<Real> M(nx_+2,0.0);
862  Real ss = dt_;
863  for (unsigned t = 0; t < nt_+1; t++) {
864  ss = ((t < nt_ && t > 0) ? dt_ : 0.5*dt_);
865  for (unsigned n = 0; n < nx_+2; n++) {
866  Z[n] = (*zp)[t*(nx_+2)+n];
867  }
868  apply_mass(M,Z);
869  for (unsigned n = 0; n < nx_+2; n++) {
870  (*gp)[t*(nx_+2)+n] = ss*alpha_*M[n];
871  }
872  }
873  }
874 
876  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
877  Teuchos::RCP<std::vector<Real> > hvp = Teuchos::rcp_const_cast<std::vector<Real> >(
878  (Teuchos::dyn_cast<const ROL::StdVector<Real> >(hv)).getVector());
879  Teuchos::RCP<const std::vector<Real> > vp =
880  (Teuchos::dyn_cast<ROL::StdVector<Real> >(const_cast<ROL::Vector<Real> &>(v))).getVector();
881  // COMPUTE GRADIENT WRT U
882  std::vector<Real> V(nx_,0.0);
883  std::vector<Real> M(nx_,0.0);
884  Real ss = 0.5*dt_;
885  for (unsigned t = 0; t < nt_; t++) {
886  ss = ((t < nt_-1) ? dt_ : 0.5*dt_);
887  for (unsigned n = 0; n < nx_; n++) {
888  V[n] = (*vp)[t*nx_+n];
889  }
890  apply_mass(M,V);
891  for (unsigned n = 0; n < nx_; n++) {
892  (*hvp)[t*nx_+n] = ss*M[n];
893  }
894  }
895  }
896 
898  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
899  hv.zero();
900  }
901 
903  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
904  hv.zero();
905  }
906 
908  const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
909  Teuchos::RCP<std::vector<Real> > hvp = Teuchos::rcp_const_cast<std::vector<Real> >(
910  (Teuchos::dyn_cast<const ROL::StdVector<Real> >(hv)).getVector());
911  Teuchos::RCP<const std::vector<Real> > vp =
912  (Teuchos::dyn_cast<ROL::StdVector<Real> >(const_cast<ROL::Vector<Real> &>(v))).getVector();
913  // COMPUTE GRADIENT WRT Z
914  std::vector<Real> V(nx_+2,0.0);
915  std::vector<Real> M(nx_+2,0.0);
916  Real ss = 0.0;
917  for (unsigned t = 0; t < nt_+1; t++) {
918  ss = ((t < nt_ && t > 0) ? dt_ : 0.5*dt_);
919  for (unsigned n = 0; n < nx_+2; n++) {
920  V[n] = (*vp)[t*(nx_+2)+n];
921  }
922  apply_mass(M,V);
923  for (unsigned n = 0; n < nx_+2; n++) {
924  (*hvp)[t*(nx_+2)+n] = ss*alpha_*M[n];
925  }
926  }
927  }
928 };
Provides the interface to evaluate simulation-based objective functions.
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:321
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:153
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 adjoint of the partial constraint Hessian at , , to vector in direction ...
Definition: example_03.hpp:726
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)
Real evaluate_target(Real x)
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 adjoint of the partial constraint Hessian at , , to vector in direction ...
Definition: example_03.hpp:730
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:572
Real dot(const std::vector< Real > &x, const std::vector< Real > &y)
Contains definitions of custom data types in ROL.
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:117
Real compute_norm(const std::vector< Real > &r)
Definition: example_03.hpp:84
EqualityConstraint_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:429
void solve(ROL::Vector< Real > &c, ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Given , solve for .
Definition: example_03.hpp:482
virtual void zero()
Set to zero vector.
Definition: ROL_Vector.hpp:157
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:74
void apply_control_jacobian(std::vector< Real > &jv, const std::vector< Real > &v, bool adjoint=false)
Definition: example_03.hpp:293
void compute_residual(std::vector< Real > &r, const std::vector< Real > &u, const std::vector< Real > &z)
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:277
Defines the equality constraint operator interface for simulation-based optimization.
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:907
Real value(const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Compute value.
Definition: example_03.hpp:806
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:833
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:608
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:178
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:633
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:516
Real dot(const std::vector< Real > &x, const std::vector< Real > &y)
Definition: example_03.hpp:762
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:902
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:897
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:854
void scale(std::vector< Real > &u, const Real alpha=0.0)
Definition: example_03.hpp:111
Objective_BurgersControl(Real alpha=1.e-4, int nx=128, int nt=100, Real T=1.0)
Definition: example_03.hpp:800
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:543
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:875
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:443
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:663
void apply_mass(std::vector< Real > &Mu, const std::vector< Real > &u)
Definition: example_03.hpp:779
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:368
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:238
virtual Real norm() const =0
Returns where .
void update(std::vector< Real > &u, const std::vector< Real > &s, const Real alpha=1.0)
Definition: example_03.hpp:105
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 adjoint of the partial constraint Hessian at , , to vector in direction ...
Definition: example_03.hpp:722
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:208
static const double ROL_EPSILON
Platform-dependent machine epsilon.
Definition: ROL_Types.hpp:118
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 adjoint of the partial constraint Hessian at , , to vector in direction ...
Definition: example_03.hpp:697