ROL
ROL_PoissonControl.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 #ifndef USE_HESSVEC
16 #define USE_HESSVEC 1
17 #endif
18 
19 #ifndef ROL_POISSONCONTROL_HPP
20 #define ROL_POISSONCONTROL_HPP
21 
22 #include "ROL_StdVector.hpp"
23 #include "ROL_TestProblem.hpp"
24 
25 namespace ROL {
26 namespace ZOO {
27 
30 template<class Real>
31 class Objective_PoissonControl : public Objective<Real> {
32 
33 typedef std::vector<Real> vector;
34 typedef Vector<Real> V;
36 
37 typedef typename vector::size_type uint;
38 
39 private:
40  Real alpha_;
41 
42  ROL::Ptr<const vector> getVector( const V& x ) {
43 
44  return dynamic_cast<const SV&>(x).getVector();
45  }
46 
47  ROL::Ptr<vector> getVector( V& x ) {
48 
49  return dynamic_cast<SV&>(x).getVector();
50  }
51 
52 public:
53 
54  Objective_PoissonControl(Real alpha = 1.e-4) : alpha_(alpha) {}
55 
56  void apply_mass(Vector<Real> &Mz, const Vector<Real> &z ) {
57 
58 
59  ROL::Ptr<const vector> zp = getVector(z);
60  ROL::Ptr<vector> Mzp = getVector(Mz);
61 
62  uint n = zp->size();
63  Real h = 1.0/((Real)n+1.0);
64  for (uint i=0; i<n; i++) {
65  if ( i == 0 ) {
66  (*Mzp)[i] = h/6.0*(4.0*(*zp)[i] + (*zp)[i+1]);
67  }
68  else if ( i == n-1 ) {
69  (*Mzp)[i] = h/6.0*((*zp)[i-1] + 4.0*(*zp)[i]);
70  }
71  else {
72  (*Mzp)[i] = h/6.0*((*zp)[i-1] + 4.0*(*zp)[i] + (*zp)[i+1]);
73  }
74  }
75  }
76 
77  void solve_poisson(Vector<Real> & u, const Vector<Real> & z) {
78 
79 
80 
81 
82  ROL::Ptr<vector> up = getVector(u);
83 
84  uint n = up->size();
85  Real h = 1.0/((Real)n+1.0);
86  SV b( ROL::makePtr<vector>(n,0.0) );
87  ROL::Ptr<vector> bp = getVector(b);
88  apply_mass(b,z);
89 
90  Real d = 2.0/h;
91  Real o = -1.0/h;
92  Real m = 0.0;
93  vector c(n,o);
94  c[0] = c[0]/d;
95  (*up)[0] = (*bp)[0]/d;
96  for ( uint i = 1; i < n; i++ ) {
97  m = 1.0/(d - o*c[i-1]);
98  c[i] = c[i]*m;
99  (*up)[i] = ( (*bp)[i] - o*(*up)[i-1] )*m;
100  }
101  for ( uint i = n-1; i > 0; i-- ) {
102  (*up)[i-1] = (*up)[i-1] - c[i-1]*(*up)[i];
103  }
104  }
105 
106  Real evaluate_target(Real x) {
107  Real val = 1.0/3.0*std::pow(x,4.0) - 2.0/3.0*std::pow(x,3.0) + 1.0/3.0*x + 8.0*alpha_;
108  return val;
109  }
110 
111  Real value( const Vector<Real> &z, Real &tol ) {
112 
113 
114 
115  ROL::Ptr<const vector> zp = getVector(z);
116  uint n = zp->size();
117  Real h = 1.0/((Real)n+1.0);
118  // SOLVE STATE EQUATION
119  SV u( ROL::makePtr<vector>(n,0.0) );
120  solve_poisson(u,z);
121  ROL::Ptr<vector> up = getVector(u);
122 
123  Real val = 0.0;
124  Real res = 0.0;
125  Real res1 = 0.0;
126  Real res2 = 0.0;
127  Real res3 = 0.0;
128  for (uint i=0; i<n; i++) {
129  res = alpha_*(*zp)[i];
130  if ( i == 0 ) {
131  res *= h/6.0*(4.0*(*zp)[i] + (*zp)[i+1]);
132  res1 = (*up)[i]-evaluate_target((Real)(i+1)*h);
133  res2 = (*up)[i+1]-evaluate_target((Real)(i+2)*h);
134  res += h/6.0*(4.0*res1 + res2)*res1;
135  }
136  else if ( i == n-1 ) {
137  res *= h/6.0*((*zp)[i-1] + 4.0*(*zp)[i]);
138  res1 = (*up)[i-1]-evaluate_target((Real)(i)*h);
139  res2 = (*up)[i]-evaluate_target((Real)(i+1)*h);
140  res += h/6.0*(res1 + 4.0*res2)*res2;
141  }
142  else {
143  res *= h/6.0*((*zp)[i-1] + 4.0*(*zp)[i] + (*zp)[i+1]);
144  res1 = (*up)[i-1]-evaluate_target((Real)(i)*h);
145  res2 = (*up)[i]-evaluate_target((Real)(i+1)*h);
146  res3 = (*up)[i+1]-evaluate_target((Real)(i+2)*h);
147  res += h/6.0*(res1 + 4.0*res2 + res3)*res2;
148  }
149  val += 0.5*res;
150  }
151  return val;
152  }
153 
154  void gradient( Vector<Real> &g, const Vector<Real> &z, Real &tol ) {
155 
156 
157 
158  ROL::Ptr<const vector> zp = getVector(z);
159  ROL::Ptr<vector> gp = getVector(g);
160 
161  uint n = zp->size();
162  Real h = 1.0/((Real)n+1.0);
163 
164  // SOLVE STATE EQUATION
165  SV u( ROL::makePtr<vector>(n,0.0) );
166  solve_poisson(u,z);
167  ROL::Ptr<vector> up = getVector(u);
168 
169  // SOLVE ADJOINT EQUATION
170  StdVector<Real> res( ROL::makePtr<std::vector<Real>>(n,0.0) );
171  ROL::Ptr<vector> rp = getVector(res);
172 
173  for (uint i=0; i<n; i++) {
174  (*rp)[i] = -((*up)[i]-evaluate_target((Real)(i+1)*h));
175  }
176 
177  SV p( ROL::makePtr<vector>(n,0.0) );
178  solve_poisson(p,res);
179  ROL::Ptr<vector> pp = getVector(p);
180 
181  Real res1 = 0.0;
182  Real res2 = 0.0;
183  Real res3 = 0.0;
184  for (uint i=0; i<n; i++) {
185  if ( i == 0 ) {
186  res1 = alpha_*(*zp)[i] - (*pp)[i];
187  res2 = alpha_*(*zp)[i+1] - (*pp)[i+1];
188  (*gp)[i] = h/6.0*(4.0*res1 + res2);
189  }
190  else if ( i == n-1 ) {
191  res1 = alpha_*(*zp)[i-1] - (*pp)[i-1];
192  res2 = alpha_*(*zp)[i] - (*pp)[i];
193  (*gp)[i] = h/6.0*(res1 + 4.0*res2);
194  }
195  else {
196  res1 = alpha_*(*zp)[i-1] - (*pp)[i-1];
197  res2 = alpha_*(*zp)[i] - (*pp)[i];
198  res3 = alpha_*(*zp)[i+1] - (*pp)[i+1];
199  (*gp)[i] = h/6.0*(res1 + 4.0*res2 + res3);
200  }
201  }
202  }
203 #if USE_HESSVEC
204  void hessVec( Vector<Real> &hv, const Vector<Real> &v, const Vector<Real> &z, Real &tol ) {
205 
206 
207 
208  ROL::Ptr<const vector> zp = getVector(z);
209  ROL::Ptr<const vector> vp = getVector(v);
210  ROL::Ptr<vector> hvp = getVector(hv);
211 
212  uint n = zp->size();
213  Real h = 1.0/((Real)n+1.0);
214 
215  // SOLVE STATE EQUATION
216  SV u( ROL::makePtr<vector>(n,0.0) );
217  solve_poisson(u,v);
218  ROL::Ptr<vector> up = getVector(u);
219 
220  // SOLVE ADJOINT EQUATION
221  SV p( ROL::makePtr<vector>(n,0.0) );
222 
223  solve_poisson(p,u);
224  ROL::Ptr<vector> pp = getVector(p);
225 
226  Real res1 = 0.0;
227  Real res2 = 0.0;
228  Real res3 = 0.0;
229  for (uint i=0; i<n; i++) {
230  if ( i == 0 ) {
231  res1 = alpha_*(*vp)[i] + (*pp)[i];
232  res2 = alpha_*(*vp)[i+1] + (*pp)[i+1];
233  (*hvp)[i] = h/6.0*(4.0*res1 + res2);
234  }
235  else if ( i == n-1 ) {
236  res1 = alpha_*(*vp)[i-1] + (*pp)[i-1];
237  res2 = alpha_*(*vp)[i] + (*pp)[i];
238  (*hvp)[i] = h/6.0*(res1 + 4.0*res2);
239  }
240  else {
241  res1 = alpha_*(*vp)[i-1] + (*pp)[i-1];
242  res2 = alpha_*(*vp)[i] + (*pp)[i];
243  res3 = alpha_*(*vp)[i+1] + (*pp)[i+1];
244  (*hvp)[i] = h/6.0*(res1 + 4.0*res2 + res3);
245  }
246  }
247  }
248 #endif
249 };
250 
251 template<class Real>
252 class getPoissonControl : public TestProblem<Real> {
253 public:
255 
256  Ptr<Objective<Real>> getObjective(void) const {
257  // Instantiate Objective Function
258  return ROL::makePtr<Objective_PoissonControl<Real>>();
259  }
260 
261  Ptr<Vector<Real>> getInitialGuess(void) const {
262  // Problem dimension
263  int n = 512;
264  // Get Initial Guess
265  ROL::Ptr<std::vector<Real> > x0p = ROL::makePtr<std::vector<Real>>(n,0.0);
266  for (int i=0; i<n; i++) {
267  (*x0p)[i] = 0.0;
268  }
269  return ROL::makePtr<StdVector<Real>>(x0p);
270  }
271 
272  Ptr<Vector<Real>> getSolution(const int i = 0) const {
273  // Problem dimension
274  int n = 512;
275  // Get Solution
276  ROL::Ptr<std::vector<Real> > xp = ROL::makePtr<std::vector<Real>>(n,0.0);
277  Real h = 1.0/((Real)n+1.0), pt = 0.0;
278  for( int i = 0; i < n; i++ ) {
279  pt = (Real)(i+1)*h;
280  (*xp)[i] = 4.0*pt*(1.0-pt);
281  }
282  return ROL::makePtr<StdVector<Real>>(xp);
283  }
284 };
285 
286 } // End ZOO Namespace
287 } // End ROL Namespace
288 
289 #endif
Provides the interface to evaluate objective functions.
typename PV< Real >::size_type size_type
Ptr< Vector< Real > > getSolution(const int i=0) const
Real value(const Vector< Real > &z, Real &tol)
Compute value.
virtual void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply Hessian approximation to vector.
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:46
Ptr< Objective< Real > > getObjective(void) const
void solve_poisson(Vector< Real > &u, const Vector< Real > &z)
void gradient(Vector< Real > &g, const Vector< Real > &z, Real &tol)
Compute gradient.
Poisson distributed control.
Contains definitions of test objective functions.
void apply_mass(Vector< Real > &Mz, const Vector< Real > &z)
ROL::Ptr< const vector > getVector(const V &x)
Ptr< Vector< Real > > getInitialGuess(void) const