ROL
ROL_Constraint_SerialSimOpt.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 
10 #pragma once
11 #ifndef ROL_CONSTRAINT_SERIALSIMOPT_HPP
12 #define ROL_CONSTRAINT_SERIALSIMOPT_HPP
13 
14 namespace ROL {
15 
40 #include "ROL_VectorWorkspace.hpp"
41 
42 
43 template<typename Real>
45 
46  using V = Vector<Real>;
49 
50  using size_type = typename PV<Real>::size_type;
51 
52 private:
53 
54  const Ptr<Con> con_; // Constraint for single time step
55  const Ptr<V> ui_; // Initial condition
56  Ptr<V> u0_; // Zero sim vector on time step
57  Pre<V> z0_; // Zero opt vector on time step
58 
59  VectorWorkspace<Real> workspace_; // Memory management
60 
61 protected:
62 
63  PV& partition( V& x ) const { return static_cast<PV&>(x); }
64 
65  const V& partition( const V& x ) const { return static_cast<const PV&>(x); }
66 
67 
68 public:
69 
70  Constraint_SerialSimOpt( const Ptr<Con>& con, const V& ui, const V& zi ) :
71  Constraint_SimOpt<Real>(), con_(con),
72  ui_( workspace_.copy( ui ) ),
73  u0_( workspace_.clone( ui ) ),
74  z0_( workspace_.clone( ui ) ) {
75  u0_->zero(); z0_->zero();
76  }
77 
78 // virtual void update_1( const V& u, const V& z, bool flag = true, int iter = -1 ) override {
79 //
80 // }
81 //
82 // virtual void update_2( const V& u, const V& z, bool flag = true, int iter = -1 ) override {
83 //
84 // }
85 
86  virtual void value( V& c, const V& u, const V& z, Real& tol ) override {
87 
88  auto& cp = partition(c);
89  auto& up = partition(u);
90  auto& zp = partition(z);
91 
92  // First interval value uses initial condition
93  con_->value( *(cp.get(0)), *ui_, *(zp.get(0)), tol );
94 
95  for( size_type k=1, k<cp.numVectors(), ++k ) {
96  con_->value( *(cp.get(k)), *(up.get(k-1)), *(up.get(k)), *(zp.get(k)), tol );
97  }
98  }
99 
100  virtual void solve( V& c, V& u, const V& z, Real& tol ) override {
102  auto& cp = partition(c);
103  auto& up = partition(u);
104  auto& zp = partition(z);
105 
106  // First interval value uses initial condition
107  con_->solve( *(cp.get(0)), *ui_, *(zp.get(0)), tol );
109  for( size_type k=1, k<cp.numVectors(), ++k ) {
110  con_->solve( *(cp.get(k)), *(up.get(k-1)), *(up.get(k)), *(zp.get(k)), tol );
111  }
112  }
113 
114  virtual void applyJacobian_1( V& jv, const V& v, const V& u, const V& z, Real& tol ) override {
115 
116  auto& jvp = partition(jv); auto& vp = partition(v);
117  auto& up = partition(u); auto& zp = partition(z);
118 
119  auto tmp = workspace_.clone( up.get(0) );
120 
121  con_->applyJacobian_1_new( *(jvp.get(0)), *(vp.get(0)), *u0_, *(up.get(0)), *(zp.get(0)), tol );
122 
123  for( size_type k=1; k<jvp.numVectors(); ++jvp ) {
124  con_->applyJacobian_1_new( *(jvp.get(k)), *(vp.get(k)), *(up.get(k)), *(up.get(k-1)), *(zp.get(k)), tol );
125  con_->applyJacobian_1_old( *tmp, *(vp.get(k)), *(up.get(k)), *(up.get(k-1)), *(zp.get(k)), tol );
126  jvp.plus(tmp);
127  }
128  } // applyJacobian_1
129 
130  virtual void applyJacobian_2( V& jv, const V& v, const V& u, const V& z, Real& tol) override {
131 
132  auto& jvp = partition(jv); auto& vp = partition(v);
133  auto& up = partition(u); auto& zp = partition(z);
134 
135  for( size_type k=0; k<jvp.numVectors(); ++k ) {
136  con_->applyJacobian_2( *(jvp.get(k)), *(vp.get(k)), *(up.get(k)), *(up.get(k-1)), *(zp.get(k)), tol );
137  }
138  } // applyJacobian_2
139 
140 
141  virtual void applyInverseJacobian_1( V& ijv, const V& v, const V& u,
142  const V& z, Real& tol ) override {
144  auto& ijvp = partition(ijv); auto& vp = partition(v);
145  auto& up = partition(u); auto& zp = partition(z);
146 
147  // First step
148  con_->applyInverseJacobian_1_new( *(ijvp.get(0)), *(vp.get(0), *u0_, *(up.get(0)), *(zp.get(0)), tol );
149 
150  auto tmp = workspace.clone( vp.get(0) );
152  for( size_type k=1; k<ijvp.numVectors(); ++k ) {
153 
154  con_->applyJacobian_1_old( *tmp, *(ijvp.get(k)), *(up.get(k-1)), *(up.get(k)), *(zp.get(k)), tol );
155  tmp->scale(-1.0);
156  tmp->plus( *(vp.get(k) );
157 
158  con_->applyInverseJacobian_1_new( *(ijvp.get(k)), *tmp, *(up.get(k-1)), *(up.get(k)), *(zp.get(k)), tol );
159  }
160  } // applyInverseJacobian_1
161 
162 
163 
164  virtual void applyAdjointJacobian_1( V& ajv, const V& v, const V& u,
165  const V& z, Real& tol) override {
166 
167  auto& ajvp = partition(ajv); auto& vp = partition(v);
168  auto& up = partition(u); auto& zp = partition(z);
169 
170  auto tmp = workspace.clone( ajvp.get(0) );
171 
172  size_type N = ajvp.numVectors()-1;
173 
174  // First step
175  con_->applyAdjointJacobian_1_new( *(ajvp.get(0)), *(vp.get(0)), *u0_, *(up.get(0)), *(zp.get(0)) );
176  con_->applyAdjointJacobian_1_old( *(ajvp.get(0)), *(vp.get(1)), *u0_, *(up.get(0)), *(zp.get(0)) );
177 
178  for( size_type k=1; k<N; ++k ) {
179  con_->applyAdjointJacobian_1_new( *(ajvp.get(k)), *(vp.get(k)), *(up.get(k-1)), *(up.get(k)), *(zp.get(k)) );
180  con_->applyAdjointJacobian_1_old( *tmp, *(vp.get(k+1)), *(up.get(k-1)), *(up.get(k)), *(zp.get(k)) );
181  ajvp.plus(*tmp);
182  }
183 
184  // Last step
185  con_->applyAdjointJacobian_1_new( *(ajvp.get(k)), *(vp.get(k)), *(up.get(k-1)), *(up.get(k)), *(zp.get(k)) );
186 
187  } // applyAdjointJacobian_1
188 
189 
190 
191  virtual void applyAdjointJacobian_1( V& ajv, const V& v, const V& u, const V &z,
192  const V& dualv, Real& tol) override {
193 
194  auto& ajvp = partition(ajv); auto& vp = partition(dualv);
195  auto& up = partition(u); auto& zp = partition(z);
196  auto tmp = workspace.clone( ajvp.get(0) );
197 
198  size_type N = ajvp.numVectors()-1;
199 
200  // First step
201  con_->applyAdjointJacobian_1_new( *(ajvp.get(0)), *(vp.get(0)), *u0_, *(up.get(0)), *(zp.get(0)) );
202  con_->applyAdjointJacobian_1_old( *(ajvp.get(0)), *(vp.get(1)), *u0_, *(up.get(0)), *(zp.get(0)) );
203 
204  for( size_type k=1; k<N; ++k ) {
205  con_->applyAdjointJacobian_1_new( *(ajvp.get(k)), *(vp.get(k)), *(up.get(k-1)), *(up.get(k)), *(zp.get(k)) );
206  con_->applyAdjointJacobian_1_old( *tmp, *(vp.get(k+1)), *(up.get(k-1)), *(up.get(k)), *(zp.get(k)) );
207  ajvp.plus(*tmp);
208  }
209 
210  // Last step
211  con_->applyAdjointJacobian_1_new( *(ajvp.get(k)), *(vp.get(k)), *(up.get(k-1)), *(up.get(k)), *(zp.get(k)) );
212 
213  } // applyAdjointJacobian_1
214 
215 
216 
217  virtual void applyAdjointJacobian_2( V& ajv, const V& v, const V& u,
218  const V& z, Real& tol) override {
219  auto& ajvp = partition(ajv); auto& vp = partition(v);
220  auto& up = partition(u); auto& zp = partition(z);
221 
222  for( size_type k=0; k<jvp.numVectors(); ++k ) {
223  con_->applyAdjointJacobian_2_time( *(jvp.get(k)), *(vp.get(k)), *(up.get(k)), *(up.get(k-1)), *(zp.get(k)), tol );
224  }
225  } // applyAdjointJacobian_2
226 
227 
228  virtual void applyAdjointJacobian_2( V& ajv, const V& v, const V& u, const V& z,
229  const V& dualv, Real& tol) override {
230 
231  auto& ajvp = partition(ajv); auto& vp = partition(v);
232  auto& up = partition(u); auto& zp = partition(z);
233 
234  for( size_type k=0; k<jvp.numVectors(); ++k ) {
235  con_->applyAdjointJacobian_2_time( *(jvp.get(k)), *(vp.get(k)), *(up.get(k)), *(up.get(k-1)), *(zp.get(k)), tol );
236  }
237  }
238 
240  virtual void applyInverseAdjointJacobian_1( V& iajv, const V& v, const V& u,
241  const V& z, Real& tol) override {
242 
243  auto& iajvp = partition(iajv); auto& vp = partition(v);
244  auto& up = partition(u); auto& zp = partition(z);
245 
246  size_type N = iajvp.numVectors()-1;
247 
248  // Last Step
249  con_->applyInverseAdjointJacobian_1_new( *(iajvp.get(N)), *(vp.get(N)), *(up.get(N-1)), *(up.get(N)), *(zp.get(N)), tol );
250 
251  auto tmp = workspace_.clone( vp.get(0) );
252 
253  for( size_type k=N-1; k>0; --k ) {
254  con_->applyAdjointJacobian_1_old( *tmp, *(iajvp.get(k+1)), *(up.get(k-1)), *(up.get(k)), *(zp.get(k)), tol );
255 
256  tmp->scale( -1.0 );
257  tmp->plus( *(vp.get(k) );
258 
259  con_->applyAdjointJacobian_1_new( *(iajvp.get(k)), *tmp, *(up.get(k-1)), *(up.get(k)), *(zp.get(k)), tol );
260  }
261 
262  // "First step"
263  con_->applyAdjointJacobian_1_old( *tmp, *(iajvp.get(1)), *u0_, *(up.get(0)), *(zp.get(0)), tol );
264 
265  tmp->scale( -1.0 );
266  tmp->plus( *(vp.get(0) );
267 
268  con_->applyAdjointJacobian_1_new( *(iajvp.get(0)), *tmp, *u0_, *(up.get(0)), *(zp.get(0)), tol );
269  }
270 
272  virtual void applyAdjointHessian_11( V& ahwv, const V& w, const V& v,
273  const V& u, const V& z, Real& tol) override {
274  // TODO
275  }
276 
277 
278  virtual void applyAdjointHessian_11( V& ahwv, const V& w, const V& v,
279  const V& u, const V& z, Real& tol) override {
280  ahwv.zero();
281  }
282 
283  virtual void applyAdjointHessian_11( V& ahwv, const V& w, const V& v,
284  const V& u, const V& z, Real& tol) override {
285  ahwv.zero();
286  }
287 
288  virtual void applyAdjointHessian_11( V& ahwv, const V& w, const V& v,
289  const V& u, const V& z, Real& tol) override {
290  ahwv.zero();
291  }
292 
293 
294 
295 
296 
297 
298 
299 
300 
301 
302 }; // class Constraint_SerialSimOpt
304 
305 
306 } // namespace ROL
307 
308 
309 #endif // ROL_CONSTRAINT_SERIALSIMOPT_HPP
310 
typename PV< Real >::size_type size_type
virtual void applyJacobian_2(V &jv, const V &v, const V &u, const V &z, Real &tol) override
Apply the partial constraint Jacobian at , , to the vector .
Defines the linear algebra of vector space on a generic partitioned vector.
virtual void applyJacobian_1(V &jv, const V &v, const V &u, const V &z, Real &tol) override
Apply the partial constraint Jacobian at , , to the vector .
virtual void applyAdjointHessian_11(V &ahwv, const V &w, const V &v, const V &u, const V &z, Real &tol) override
Apply the simulation-space derivative of the adjoint of the constraint simulation-space Jacobian at ...
virtual void zero()
Set to zero vector.
Definition: ROL_Vector.hpp:133
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:46
Defines the time dependent constraint operator interface for simulation-based optimization.
virtual void applyAdjointJacobian_1(V &ajv, const V &v, const V &u, const V &z, Real &tol) override
Apply the adjoint of the partial constraint Jacobian at , , to the vector . This is the primary inter...
virtual void applyAdjointJacobian_2(V &ajv, const V &v, const V &u, const V &z, Real &tol) override
Apply the adjoint of the partial constraint Jacobian at , , to vector . This is the primary interface...
Unifies the constraint defined on a single time step that are defined through the Constraint_TimeSimO...
const V & partition(const V &x) const
virtual void applyAdjointJacobian_2(V &ajv, const V &v, const V &u, const V &z, const V &dualv, Real &tol) override
Apply the adjoint of the partial constraint Jacobian at , , to vector . This is the secondary interfa...
typename PV< Real >::size_type size_type
std::vector< PV >::size_type size_type
virtual void applyInverseAdjointJacobian_1(V &iajv, const V &v, const V &u, const V &z, Real &tol) override
Apply the inverse of the adjoint of the partial constraint Jacobian at , , to the vector ...
virtual void applyInverseJacobian_1(V &ijv, const V &v, const V &u, const V &z, Real &tol) override
Apply the inverse partial constraint Jacobian at , , to the vector .
virtual void solve(V &c, V &u, const V &z, Real &tol) override
Given , solve for .
Constraint_SerialSimOpt(const Ptr< Con > &con, const V &ui, const V &zi)
virtual void value(V &c, const V &u, const V &z, Real &tol) override
Evaluate the constraint operator at .