ROL
gross-pitaevskii/example_01.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
70#include <iostream>
71
72#include "ROL_Stream.hpp"
73#include "Teuchos_GlobalMPISession.hpp"
74
75#include "ROL_ParameterList.hpp"
76#include "ROL_StdVector.hpp"
77#include "ROL_Objective.hpp"
78#include "ROL_Constraint.hpp"
79#include "ROL_Algorithm.hpp"
80#include "ROL_CompositeStep.hpp"
82
83
84using namespace ROL;
85
86template<class Real>
88
89 typedef std::vector<Real> vector;
90 typedef Vector<Real> V;
92
93 typedef typename vector::size_type uint;
94
95private:
96
97
99 Real g_;
100
103
105 Real dx_;
106
108 ROL::Ptr<const vector> Vp_;
109
110 ROL::Ptr<const vector> getVector( const V& x ) {
111
112 return dynamic_cast<const SV&>(x).getVector();
113 }
114
115 ROL::Ptr<vector> getVector( V& x ) {
116
117 return dynamic_cast<SV&>(x).getVector();
118 }
119
120
121
123
125 void applyK(const Vector<Real> &v, Vector<Real> &kv) {
126
127 using namespace Teuchos;
128
129 // Pointer to direction vector
130 ROL::Ptr<const vector> vp = getVector(v);
131
132 // Pointer to action of Hessian on direction vector
133 ROL::Ptr<vector> kvp = getVector(kv);
134
135 Real dx2 = dx_*dx_;
136
137 (*kvp)[0] = (2.0*(*vp)[0]-(*vp)[1])/dx2;
138
139 for(uint i=1;i<nx_-1;++i) {
140 (*kvp)[i] = (2.0*(*vp)[i]-(*vp)[i-1]-(*vp)[i+1])/dx2;
141 }
142
143 (*kvp)[nx_-1] = (2.0*(*vp)[nx_-1]-(*vp)[nx_-2])/dx2;
144 }
145
146 public:
147
148 Objective_GrossPitaevskii(const Real &g, const Vector<Real> &V) : g_(g),
149 Vp_(getVector(V)) {
150 nx_ = Vp_->size();
151 dx_ = (1.0/(1.0+nx_));
152 }
153
155
159 Real value( const Vector<Real> &psi, Real &tol ) {
160
161
162
163 // Pointer to opt vector
164 ROL::Ptr<const vector> psip = getVector(psi);
165
166 // Pointer to K applied to opt vector
167 ROL::Ptr<V> kpsi = psi.clone();
168 ROL::Ptr<vector> kpsip = getVector(*kpsi);
169
170 Real J = 0;
171
172 this->applyK(psi,*kpsi);
173
174 for(uint i=0;i<nx_;++i) {
175 J += (*psip)[i]*(*kpsip)[i] + (*Vp_)[i]*pow((*psip)[i],2) + g_*pow((*psip)[i],4);
176 }
177
178 // Rescale for numerical integration by trapezoidal rule
179 J *= 0.5*dx_;
180
181 return J;
182 }
183
185
186 void gradient( Vector<Real> &g, const Vector<Real> &psi, Real &tol ) {
187
188
189
190 // Pointer to opt vector
191 ROL::Ptr<const vector> psip = getVector(psi);
192
193 // Pointer to gradient vector
194 ROL::Ptr<vector> gp = getVector(g);
195
196 // Pointer to K applied to opt vector
197 ROL::Ptr<V> kpsi = psi.clone();
198 ROL::Ptr<vector> kpsip = getVector(*kpsi);
199
200 applyK(psi,*kpsi);
201
202 for(uint i=0;i<nx_;++i) {
203 (*gp)[i] = ((*kpsip)[i] + (*Vp_)[i]*(*psip)[i] + 2.0*g_*pow((*psip)[i],3))*dx_;
204 }
205
206 }
207
208
209
211
212 void hessVec( Vector<Real> &hv, const Vector<Real> &v, const Vector<Real> &psi, Real &tol ) {
213
214
215
216 // Pointer to opt vector
217 ROL::Ptr<const vector> psip = getVector(psi);
218
219 // Pointer to direction vector
220 ROL::Ptr<const vector> vp = getVector(v);
221
222 // Pointer to action of Hessian on direction vector
223 ROL::Ptr<vector> hvp = getVector(hv);
224
225 applyK(v,hv);
226
227 for(uint i=0;i<nx_;++i) {
228 (*hvp)[i] *= dx_;
229 (*hvp)[i] += ( (*Vp_)[i] + 6.0*g_*pow((*psip)[i],2) )*(*vp)[i]*dx_;
230 }
231 }
232
233};
234
235
236template<class Real>
238
239 typedef std::vector<Real> vector;
242
243 typedef typename vector::size_type uint;
244
245
246private:
248 Real dx_;
249
250 ROL::Ptr<const vector> getVector( const V& x ) {
251
252 return dynamic_cast<const SV&>(x).getVector();
253 }
254
255 ROL::Ptr<vector> getVector( V& x ) {
256
257 return dynamic_cast<SV&>(x).getVector();
258 }
259
260public:
261 Normalization_Constraint(int n, Real dx) : nx_(n), dx_(dx) {}
262
264
268 void value(Vector<Real> &c, const Vector<Real> &psi, Real &tol){
269
270
271
272 // Pointer to constraint vector (only one element)
273 ROL::Ptr<vector> cp = getVector(c);
274
275 // Pointer to optimization vector
276 ROL::Ptr<const vector> psip = getVector(psi);
277
278 (*cp)[0] = -1.0;
279 for(uint i=0;i<nx_;++i) {
280 (*cp)[0] += dx_*pow((*psip)[i],2);
281 }
282 }
283
285
287 void applyJacobian(Vector<Real> &jv, const Vector<Real> &v, const Vector<Real> &psi, Real &tol){
288
289
290
291 // Pointer to action of Jacobian of constraint on direction vector (yields scalar)
292 ROL::Ptr<vector> jvp = getVector(jv);
293
294 // Pointer to direction vector
295 ROL::Ptr<const vector> vp = getVector(v);
296
297 // Pointer to optimization vector
298 ROL::Ptr<const vector> psip = getVector(psi);
299
300 (*jvp)[0] = 0;
301 for(uint i=0;i<nx_;++i) {
302 (*jvp)[0] += 2.0*dx_*(*psip)[i]*(*vp)[i];
303 }
304 }
305
307
309 void applyAdjointJacobian(Vector<Real> &ajv, const Vector<Real> &v, const Vector<Real> &psi, Real &tol){
310
311
312
313 // Pointer to action of adjoint of Jacobian of constraint on direction vector (yields vector)
314 ROL::Ptr<vector> ajvp = getVector(ajv);
315
316 // Pointer to direction vector
317 ROL::Ptr<const vector> vp = getVector(v);
318
319 // Pointer to optimization vector
320 ROL::Ptr<const vector> psip = getVector(psi);
321
322 for(uint i=0;i<nx_;++i) {
323 (*ajvp)[i] = 2.0*dx_*(*psip)[i]*(*vp)[0];
324 }
325 }
326
328
332 const Vector<Real> &psi, Real &tol){
333
334
335
336 // The pointer to action of constraint Hessian in u,v inner product
337 ROL::Ptr<vector> ahuvp = getVector(ahuv);
338
339 // Pointer to direction vector u
340 ROL::Ptr<const vector> up = getVector(u);
341
342 // Pointer to direction vector v
343 ROL::Ptr<const vector> vp = getVector(v);
344
345 // Pointer to optimization vector
346 ROL::Ptr<const vector> psip = getVector(psi);
347
348 for(uint i=0;i<nx_;++i) {
349 (*ahuvp)[i] = 2.0*dx_*(*vp)[i]*(*up)[0];
350 }
351 }
352};
353
354
Defines a no-output stream class ROL::NullStream and a function makeStreamPtr which either wraps a re...
void applyJacobian(Vector< Real > &jv, const Vector< Real > &v, const Vector< Real > &psi, Real &tol)
Evaluate .
void value(Vector< Real > &c, const Vector< Real > &psi, Real &tol)
Evaluate .
void applyAdjointHessian(Vector< Real > &ahuv, const Vector< Real > &u, const Vector< Real > &v, const Vector< Real > &psi, Real &tol)
Evaluate .
ROL::Ptr< const vector > getVector(const V &x)
void applyAdjointJacobian(Vector< Real > &ajv, const Vector< Real > &v, const Vector< Real > &psi, Real &tol)
Evaluate .
void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &psi, Real &tol)
Evaluate .
Objective_GrossPitaevskii(const Real &g, const Vector< Real > &V)
ROL::Ptr< const vector > getVector(const V &x)
void applyK(const Vector< Real > &v, Vector< Real > &kv)
Apply finite difference operator.
void gradient(Vector< Real > &g, const Vector< Real > &psi, Real &tol)
Evaluate .
Real value(const Vector< Real > &psi, Real &tol)
Evaluate .
Defines the general constraint operator interface.
Provides the interface to evaluate objective functions.
Provides the ROL::Vector interface for scalar values, to be used, for example, with scalar constraint...
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:84
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.