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_TestProblem.hpp"
58
59namespace ROL {
60namespace ZOO {
61
64template<class Real>
65class Objective_PoissonControl : public Objective<Real> {
66
67typedef std::vector<Real> vector;
70
71typedef typename vector::size_type uint;
72
73private:
74 Real alpha_;
75
76 ROL::Ptr<const vector> getVector( const V& x ) {
77
78 return dynamic_cast<const SV&>(x).getVector();
79 }
80
81 ROL::Ptr<vector> getVector( V& x ) {
82
83 return dynamic_cast<SV&>(x).getVector();
84 }
85
86public:
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
93 ROL::Ptr<const vector> zp = getVector(z);
94 ROL::Ptr<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
112
113
114
115
116 ROL::Ptr<vector> up = getVector(u);
117
118 uint n = up->size();
119 Real h = 1.0/((Real)n+1.0);
120 SV b( ROL::makePtr<vector>(n,0.0) );
121 ROL::Ptr<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
148
149 ROL::Ptr<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( ROL::makePtr<vector>(n,0.0) );
154 solve_poisson(u,z);
155 ROL::Ptr<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
191
192 ROL::Ptr<const vector> zp = getVector(z);
193 ROL::Ptr<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( ROL::makePtr<vector>(n,0.0) );
200 solve_poisson(u,z);
201 ROL::Ptr<vector> up = getVector(u);
202
203 // SOLVE ADJOINT EQUATION
204 StdVector<Real> res( ROL::makePtr<std::vector<Real>>(n,0.0) );
205 ROL::Ptr<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( ROL::makePtr<vector>(n,0.0) );
212 solve_poisson(p,res);
213 ROL::Ptr<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
241
242 ROL::Ptr<const vector> zp = getVector(z);
243 ROL::Ptr<const vector> vp = getVector(v);
244 ROL::Ptr<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( ROL::makePtr<vector>(n,0.0) );
251 solve_poisson(u,v);
252 ROL::Ptr<vector> up = getVector(u);
253
254 // SOLVE ADJOINT EQUATION
255 SV p( ROL::makePtr<vector>(n,0.0) );
256
257 solve_poisson(p,u);
258 ROL::Ptr<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
285template<class Real>
286class getPoissonControl : public TestProblem<Real> {
287public:
289
290 Ptr<Objective<Real>> getObjective(void) const {
291 // Instantiate Objective Function
292 return ROL::makePtr<Objective_PoissonControl<Real>>();
293 }
294
295 Ptr<Vector<Real>> getInitialGuess(void) const {
296 // Problem dimension
297 int n = 512;
298 // Get Initial Guess
299 ROL::Ptr<std::vector<Real> > x0p = ROL::makePtr<std::vector<Real>>(n,0.0);
300 for (int i=0; i<n; i++) {
301 (*x0p)[i] = 0.0;
302 }
303 return ROL::makePtr<StdVector<Real>>(x0p);
304 }
305
306 Ptr<Vector<Real>> getSolution(const int i = 0) const {
307 // Problem dimension
308 int n = 512;
309 // Get Solution
310 ROL::Ptr<std::vector<Real> > xp = ROL::makePtr<std::vector<Real>>(n,0.0);
311 Real h = 1.0/((Real)n+1.0), pt = 0.0;
312 for( int i = 0; i < n; i++ ) {
313 pt = (Real)(i+1)*h;
314 (*xp)[i] = 4.0*pt*(1.0-pt);
315 }
316 return ROL::makePtr<StdVector<Real>>(xp);
317 }
318};
319
320} // End ZOO Namespace
321} // End ROL Namespace
322
323#endif
Contains definitions of test objective functions.
Provides the interface to evaluate objective functions.
virtual void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply Hessian approximation to vector.
Provides the ROL::Vector interface for scalar values, to be used, for example, with scalar constraint...
Defines the linear algebra or vector space interface.
ROL::Ptr< const vector > getVector(const V &x)
void apply_mass(Vector< Real > &Mz, const Vector< Real > &z)
void gradient(Vector< Real > &g, const Vector< Real > &z, Real &tol)
Compute gradient.
Real value(const Vector< Real > &z, Real &tol)
Compute value.
void solve_poisson(Vector< Real > &u, const Vector< Real > &z)
Ptr< Objective< Real > > getObjective(void) const
Ptr< Vector< Real > > getInitialGuess(void) const
Ptr< Vector< Real > > getSolution(const int i=0) const