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