ROL
poisson-control/example_02.hpp
Go to the documentation of this file.
1#include <vector>
2#include <cmath>
3#include <limits>
4
5#include "Teuchos_LAPACK.hpp"
6#include "Teuchos_SerialDenseMatrix.hpp"
7
8#include "ROL_Types.hpp"
9#include "ROL_StdVector.hpp"
12
13
14template<class Real>
15class FEM {
16private:
17 int nx_; // Number of intervals
18 Real kl_; // Left diffusivity
19 Real kr_; // Right diffusivity
20
21 std::vector<Real> pts_;
22 std::vector<Real> wts_;
23public:
24
25 FEM(int nx = 128, Real kl = 0.1, Real kr = 10.0) : nx_(nx), kl_(kl), kr_(kr) {
26 // Set quadrature points
27 pts_.clear();
28 pts_.resize(5,0.0);
29 pts_[0] = -0.9061798459386640;
30 pts_[1] = -0.5384693101056831;
31 pts_[2] = 0.0;
32 pts_[3] = 0.5384693101056831;
33 pts_[4] = 0.9061798459386640;
34 wts_.clear();
35 wts_.resize(5,0.0);
36 wts_[0] = 0.2369268850561891;
37 wts_[1] = 0.4786286704993665;
38 wts_[2] = 0.5688888888888889;
39 wts_[3] = 0.4786286704993665;
40 wts_[4] = 0.2369268850561891;
41 // Scale and normalize points and weights
42 Real sum = 0.0;
43 for (int i = 0; i < 5; i++) {
44 pts_[i] = 0.5*pts_[i] + 0.5;
45 sum += wts_[i];
46 }
47 for (int i = 0; i < 5; i++) {
48 wts_[i] /= sum;
49 }
50 }
51
52 int nu(void) { return nx_-1; }
53 int nz(void) { return nx_+1; }
54
55 void build_mesh(std::vector<Real> &x, const std::vector<Real> &param) {
56 // Partition mesh:
57 // nx1 is the # of intervals on the left of xp -- n1 is the # of interior nodes
58 // nx2 is the # of intervals on the right of xp -- n2 is the # of interior nodes
59 Real xp = 0.1*param[0];
60 int frac = nx_/2;
61 int rem = nx_%2;
62 int nx1 = frac;
63 int nx2 = frac;
64 if ( rem == 1 ) {
65 if (xp > 0.0) {
66 nx2++;
67 }
68 else {
69 nx1++;
70 }
71 }
72 int n1 = nx1-1;
73 int n2 = nx2-1;
74 // Compute mesh spacing
75 Real dx1 = (xp+1.0)/(Real)nx1;
76 Real dx2 = (1.0-xp)/(Real)nx2;
77 // Build uniform meshes on [-1,xp] u [xp,1]
78 x.clear();
79 x.resize(n1+n2+3,0.0);
80 for (int k = 0; k < n1+n2+3; k++) {
81 x[k] = ((k < nx1) ? (Real)k*dx1-1.0 : ((k==nx1) ? xp : (Real)(k-nx1)*dx2+xp));
82 //std::cout << x[k] << "\n";
83 }
84 }
85
86 void build_mesh(std::vector<Real> &x) {
87 int frac = nz()/16;
88 int rem = nz()%16;
89 int nz1 = frac*4;
90 int nz2 = frac;
91 int nz3 = frac*9;
92 int nz4 = frac*2;
93 for ( int i = 0; i < rem; i++ ) {
94 if ( i%4 == 0 ) { nz3++; }
95 else if ( i%4 == 1 ) { nz1++; }
96 else if ( i%4 == 2 ) { nz4++; }
97 else if ( i%4 == 3 ) { nz2++; }
98 }
99 x.clear();
100 x.resize(nz(),0.0);
101 for (int k = 0; k < nz(); k++) {
102 if ( k < nz1 ) {
103 x[k] = 0.25*(Real)k/(Real)(nz1-1) - 1.0;
104 }
105 if ( k >= nz1 && k < nz1+nz2 ) {
106 x[k] = 0.5*(Real)(k-nz1+1)/(Real)nz2 - 0.75;
107 }
108 if ( k >= nz1+nz2 && k < nz1+nz2+nz3 ) {
109 x[k] = 0.5*(Real)(k-nz1-nz2+1)/(Real)nz3 - 0.25;
110 }
111 if ( k >= nz1+nz2+nz3-1 ) {
112 x[k] = 0.75*(Real)(k-nz1-nz2-nz3+1)/(Real)nz4 + 0.25;
113 }
114 //std::cout << x[k] << "\n";
115 }
116 //x.clear();
117 //x.resize(nz(),0.0);
118 //for (int k = 0; k < nz(); k++) {
119 // x[k] = 2.0*(Real)k/(Real)(nz()-1)-1.0;
120 // //std::cout << xz[k] << "\n";
121 //}
122 }
123
124 // Build force.
125 void build_force(std::vector<Real> &F, const std::vector<Real> &param) {
126 // Build mesh
127 std::vector<Real> x(nu()+2,0.0);
128 build_mesh(x,param);
129 // Build force term
130 F.clear();
131 F.resize(x.size()-2,0.0);
132 Real pt = 0.0;
133 int size = static_cast<int>(x.size());
134 for (int i = 0; i < size-2; i++) {
135 for (int j = 0; j < 5; j++) {
136 // Integrate over [xl,x0]
137 pt = (x[i+1]-x[i])*pts_[j] + x[i];
138 F[i] += wts_[j]*(pt-x[i])*std::exp(-std::pow(pt-0.5*param[1],2.0));
139 // Integrate over [x0,xu]
140 pt = (x[i+2]-x[i+1])*pts_[j] + x[i+1];
141 F[i] += wts_[j]*(x[i+2]-pt)*std::exp(-std::pow(pt-0.5*param[1],2.0));
142 }
143 }
144 }
145
146 // Build the PDE residual Jacobian.
147 void build_jacobian_1(std::vector<Real> &dl, std::vector<Real> &d, std::vector<Real> &du,
148 const std::vector<Real> &param) {
149 // Build mesh
150 std::vector<Real> x(nu()+2,0.0);
151 build_mesh(x,param);
152 // Fill diagonal
153 int xsize = static_cast<int>(x.size());
154 d.clear();
155 d.resize(xsize-2,0.0);
156 for ( int i = 0; i < xsize-2; i++ ) {
157 if ( x[i+1] < 0.1*param[0] ) {
158 d[i] = kl_/(x[i+1]-x[i]) + kl_/(x[i+2]-x[i+1]);
159 }
160 else if ( x[i+1] > 0.1*param[0] ) {
161 d[i] = kr_/(x[i+1]-x[i]) + kr_/(x[i+2]-x[i+1]);
162 }
163 else {
164 d[i] = kl_/(x[i+1]-x[i]) + kr_/(x[i+2]-x[i+1]);
165 }
166 //std::cout << d[i] << "\n";
167 }
168 // Fill off-diagonal
169 dl.clear();
170 dl.resize(xsize-3,0.0);
171 for ( int i = 0; i < xsize-3; i++ ) {
172 if ( x[i+2] <= 0.1*param[0] ) {
173 dl[i] = -kl_/(x[i+2]-x[i+1]);
174 }
175 else if ( x[i+2] > 0.1*param[0] ) {
176 dl[i] = -kr_/(x[i+2]-x[i+1]);
177 }
178 //std::cout << dl[i] << "\n";
179 }
180 du.clear();
181 du.assign(dl.begin(),dl.end());
182 }
183
184 // Apply the PDE residual Jacobian.
185 void apply_jacobian_1(std::vector<Real> &jv, const std::vector<Real> &v, const std::vector<Real> &param) {
186 // Build Jacobian
187 std::vector<Real> dl(nu()-1,0.0);
188 std::vector<Real> d(nu(),0.0);
189 std::vector<Real> du(nu()-1,0.0);
190 build_jacobian_1(dl,d,du,param);
191 // Apply Jacobian
192 jv.clear();
193 jv.resize(d.size(),0.0);
194 int size = static_cast<int>(d.size());
195 for ( int i = 0; i < size; ++i ) {
196 jv[i] = d[i]*v[i];
197 jv[i] += ((i>0) ? dl[i-1]*v[i-1] : 0.0);
198 jv[i] += ((i<size-1) ? du[i]*v[i+1] : 0.0);
199 }
200 }
201
202 void apply_jacobian_2(std::vector<Real> &jv, const std::vector<Real> &v, const std::vector<Real> &param) {
203 // Build u mesh
204 std::vector<Real> xu(nu()+2,0.0);
205 build_mesh(xu,param);
206 // Build z mesh
207 std::vector<Real> xz(nz(),0.0);
208 build_mesh(xz);
209 // Build union of u and z meshes
210 std::vector<Real> x(xu.size()+xz.size(),0.0);
211 typename std::vector<Real>::iterator it;
212 it = std::set_union(xu.begin(),xu.end(),xz.begin(),xz.end(),x.begin());
213 x.resize(it-x.begin());
214 // Interpolate v onto x basis
215 std::vector<Real> vi;
216 int xzsize = static_cast<int>(xz.size());
217 for (it=x.begin(); it!=x.end(); it++) {
218 for (int i = 0; i < xzsize-1; i++) {
219 if ( *it >= xz[i] && *it <= xz[i+1] ) {
220 vi.push_back(v[i+1]*(*it-xz[i])/(xz[i+1]-xz[i]) + v[i]*(xz[i+1]-*it)/(xz[i+1]-xz[i]));
221 break;
222 }
223 }
224 }
225 int xsize = static_cast<int>(x.size());
226 // Apply x basis mass matrix to interpolated v
227 std::vector<Real> Mv(xsize,0.0);
228 for (int i = 0; i < xsize; i++) {
229 Mv[i] += ((i>0) ? (x[i]-x[i-1])/6.0*vi[i-1] + (x[i]-x[i-1])/3.0*vi[i] : 0.0);
230 Mv[i] += ((i<xsize-1) ? (x[i+1]-x[i])/3.0*vi[i] + (x[i+1]-x[i])/6.0*vi[i+1] : 0.0);
231 }
232 // Reduced mass times v to u basis
233 int xusize = static_cast<int>(xu.size());
234 jv.clear();
235 jv.resize(xusize-2,0.0);
236 for (int i = 0; i < xusize-2; i++) {
237 for (int j = 0; j < xsize-1; j++) {
238 jv[i] += (((x[j]>=xu[i ]) && (x[j]< xu[i+1])) ? Mv[j]*(x[j]-xu[i ])/(xu[i+1]-xu[i ]) : 0.0);
239 jv[i] += (((x[j]>=xu[i+1]) && (x[j]<=xu[i+2])) ? Mv[j]*(xu[i+2]-x[j])/(xu[i+2]-xu[i+1]) : 0.0);
240 if ( x[j] > xu[i+2] ) { break; }
241 }
242 }
243 }
244
245 void apply_adjoint_jacobian_2(std::vector<Real> &jv, const std::vector<Real> &v,
246 const std::vector<Real> &param){
247 // Build u mesh
248 std::vector<Real> xu(nu()+2,0.0);
249 build_mesh(xu,param);
250 // Build z mesh
251 std::vector<Real> xz(nz(),0.0);
252 build_mesh(xz);
253 // Build union of u and z meshes
254 std::vector<Real> x(xu.size()+xz.size(),0.0);
255 typename std::vector<Real>::iterator it;
256 it = std::set_union(xu.begin(),xu.end(),xz.begin(),xz.end(),x.begin());
257 x.resize(it-x.begin());
258 // Interpolate v onto x basis
259 std::vector<Real> vi;
260 Real val = 0.0;
261 int xusize = static_cast<int>(xu.size());
262 for (it=x.begin(); it!=x.end(); it++) {
263 for (int i = 0; i < xusize-1; i++) {
264 if ( *it >= xu[i] && *it <= xu[i+1] ) {
265 val = 0.0;
266 val += ((i!=0 ) ? v[i-1]*(xu[i+1]-*it)/(xu[i+1]-xu[i]) : 0.0);
267 val += ((i!=xusize-2) ? v[i ]*(*it-xu[i ])/(xu[i+1]-xu[i]) : 0.0);
268 vi.push_back(val);
269 break;
270 }
271 }
272 }
273 // Apply x basis mass matrix to interpolated v
274 int xsize = static_cast<int>(x.size());
275 std::vector<Real> Mv(xsize,0.0);
276 for (int i = 0; i < xsize; i++) {
277 Mv[i] += ((i>0) ? (x[i]-x[i-1])/6.0*vi[i-1] + (x[i]-x[i-1])/3.0*vi[i] : 0.0);
278 Mv[i] += ((i<xsize-1) ? (x[i+1]-x[i])/3.0*vi[i] + (x[i+1]-x[i])/6.0*vi[i+1] : 0.0);
279 }
280 // Reduced mass times v to u basis
281 jv.clear();
282 jv.resize(nz(),0.0);
283 int xzsize = static_cast<int>(xz.size());
284 for (int i = 0; i < xzsize; i++) {
285 for (int j = 0; j < xsize; j++) {
286 if ( i==0 ) {
287 jv[i] += (((x[j]>=xz[i ]) && (x[j]<=xz[i+1])) ? Mv[j]*(xz[i+1]-x[j])/(xz[i+1]-xz[i ]) : 0.0);
288 if ( x[j] > xz[i+1] ) { break; }
289 }
290 else if ( i == xzsize-1 ) {
291 jv[i] += (((x[j]>=xz[i-1]) && (x[j]<=xz[i ])) ? Mv[j]*(x[j]-xz[i-1])/(xz[i ]-xz[i-1]) : 0.0);
292 }
293 else {
294 jv[i] += (((x[j]>=xz[i-1]) && (x[j]< xz[i ])) ? Mv[j]*(x[j]-xz[i-1])/(xz[i ]-xz[i-1]) : 0.0);
295 jv[i] += (((x[j]>=xz[i ]) && (x[j]<=xz[i+1])) ? Mv[j]*(xz[i+1]-x[j])/(xz[i+1]-xz[i ]) : 0.0);
296 if ( x[j] > xz[i+1] ) { break; }
297 }
298 }
299 }
300 }
301
302 void apply_mass_state(std::vector<Real> &Mv, const std::vector<Real> &v, const std::vector<Real> &param) {
303 // Build u mesh
304 std::vector<Real> x;
305 build_mesh(x,param);
306 // Apply mass matrix
307 int size = static_cast<int>(x.size());
308 Mv.clear();
309 Mv.resize(size-2,0.0);
310 for (int i = 0; i < size-2; i++) {
311 Mv[i] = (x[i+2]-x[i])/3.0*v[i];
312 if ( i > 0 ) {
313 Mv[i] += (x[i+1]-x[i])/6.0*v[i-1];
314 }
315 if ( i < size-3 ) {
316 Mv[i] += (x[i+2]-x[i+1])/6.0*v[i+1];
317 }
318 }
319 }
320
321 void apply_mass_control(std::vector<Real> &Mv, const std::vector<Real> &v) {
322 // Build z mesh
323 std::vector<Real> x(nz(),0.0);
324 build_mesh(x);
325 // Apply mass matrix
326 int xsize = static_cast<Real>(x.size());
327 Mv.clear();
328 Mv.resize(xsize,0.0);
329 for (int i = 0; i < xsize; i++) {
330 if ( i > 0 ) {
331 Mv[i] += (x[i]-x[i-1])/6.0*v[i-1] + (x[i]-x[i-1])/3.0*v[i];
332 }
333 if ( i < xsize-1 ) {
334 Mv[i] += (x[i+1]-x[i])/3.0*v[i] + (x[i+1]-x[i])/6.0*v[i+1];
335 }
336 }
337 }
338
339 void apply_inverse_mass_control(std::vector<Real> &Mv, const std::vector<Real> &v) {
340 // Build z mesh
341 std::vector<Real> x(nz(),0.0);
342 build_mesh(x);
343 // Build mass matrix
344 std::vector<Real> d(nz(),0.0);
345 std::vector<Real> dl(nz()-1,0.0);
346 std::vector<Real> du(nz()-1,0.0);
347 for (int i = 0; i < nz(); i++) {
348 if ( i > 0 ) {
349 dl[i-1] = (x[i]-x[i-1])/6.0;
350 d[i] += (x[i]-x[i-1])/3.0;
351 }
352 if ( i < nz()-1 ) {
353 d[i] += (x[i+1]-x[i])/3.0;
354 du[i] = (x[i+1]-x[i])/6.0;
355 }
356 }
357 // Solve linear system
358 linear_solve(Mv,dl,d,du,v,false);
359 }
360
361 // Solve linear systems in which the matrix is triadiagonal.
362 // This function uses LAPACK routines for:
363 // 1.) Computing the LU factorization of a tridiagonal matrix.
364 // 2.) Use LU factors to solve the linear system.
365 void linear_solve(std::vector<Real> &u, std::vector<Real> &dl, std::vector<Real> &d,
366 std::vector<Real> &du, const std::vector<Real> &r, const bool transpose = false) {
367 u.clear();
368 u.assign(r.begin(),r.end());
369 // Perform LDL factorization
370 Teuchos::LAPACK<int,Real> lp;
371 std::vector<Real> du2(r.size()-2,0.0);
372 std::vector<int> ipiv(r.size(),0);
373 int n = r.size();
374 int info;
375 int ldb = n;
376 int nhrs = 1;
377 lp.GTTRF(n,&dl[0],&d[0],&du[0],&du2[0],&ipiv[0],&info);
378 char trans = 'N';
379 if ( transpose ) {
380 trans = 'T';
381 }
382 lp.GTTRS(trans,n,nhrs,&dl[0],&d[0],&du[0],&du2[0],&ipiv[0],&u[0],ldb,&info);
383 }
384
385 // Evaluates the difference between the solution to the PDE.
386 // and the desired profile
387 Real evaluate_target(int i, const std::vector<Real> &param) {
388 std::vector<Real> x;
389 build_mesh(x,param);
390 Real val = 0.0;
391 int example = 2;
392 switch (example) {
393 case 1: val = ((x[i]<0.5) ? 1.0 : 0.0); break;
394 case 2: val = 1.0; break;
395 case 3: val = std::abs(std::sin(8.0*M_PI*x[i])); break;
396 case 4: val = std::exp(-0.5*(x[i]-0.5)*(x[i]-0.5)); break;
397 }
398 return val;
399 }
400};
401
402template<class Real>
404private:
405 const ROL::Ptr<FEM<Real> > FEM_;
407
408/***************************************************************/
409/********** BEGIN PRIVATE MEMBER FUNCTION DECLARATION **********/
410/***************************************************************/
411 // Returns u <- (u + alpha*s) where u, s are vectors and
412 // alpha is a scalar.
413 void plus(std::vector<Real> &u, const std::vector<Real> &s, const Real alpha=1.0) {
414 int size = static_cast<int>(u.size());
415 for (int i=0; i<size; i++) {
416 u[i] += alpha*s[i];
417 }
418 }
419
420 // Returns u <- alpha*u where u is a vector and alpha is a scalar
421 void scale(std::vector<Real> &u, const Real alpha=0.0) {
422 int size = static_cast<int>(u.size());
423 for (int i=0; i<size; i++) {
424 u[i] *= alpha;
425 }
426 }
427/*************************************************************/
428/********** END PRIVATE MEMBER FUNCTION DECLARATION **********/
429/*************************************************************/
430
431public:
433
434 int getNumSolves(void) const {
435 return num_solves_;
436 }
437
438 void value(ROL::Vector<Real> &c, const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
439 ROL::Ptr<std::vector<Real> > cp = dynamic_cast<ROL::StdVector<Real>&>(c).getVector();
440 ROL::Ptr<const std::vector<Real> > up = dynamic_cast<const ROL::StdVector<Real>&>(u).getVector();
441 ROL::Ptr<const std::vector<Real> > zp = dynamic_cast<const ROL::StdVector<Real>&>(z).getVector();
442 // Apply PDE Jacobian
443 FEM_->apply_jacobian_1(*cp, *up, this->getParameter());
444 // Get Right Hand Side
445 std::vector<Real> F(FEM_->nu(),0.0);
446 FEM_->build_force(F,this->getParameter());
447 std::vector<Real> Bz(FEM_->nu(),0.0);
448 FEM_->apply_jacobian_2(Bz,*zp,this->getParameter());
449 plus(F,Bz);
450 // Add Right Hand Side to PDE Term
451 plus(*cp,F,-1.0);
452 }
453
454 void solve(ROL::Vector<Real> &c, ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
455 ROL::Ptr<std::vector<Real> > cp = dynamic_cast<ROL::StdVector<Real>&>(c).getVector();
456 ROL::Ptr<std::vector<Real> > up = dynamic_cast<ROL::StdVector<Real>&>(u).getVector();
457 ROL::Ptr<const std::vector<Real> > zp = dynamic_cast<const ROL::StdVector<Real>&>(z).getVector();
458 // Get PDE Jacobian
459 std::vector<Real> dl(FEM_->nu()-1,0.0);
460 std::vector<Real> d(FEM_->nu(),0.0);
461 std::vector<Real> du(FEM_->nu()-1,0.0);
462 FEM_->build_jacobian_1(dl,d,du,this->getParameter());
463 // Get Right Hand Side
464 std::vector<Real> F(FEM_->nu(),0.0);
465 FEM_->build_force(F,this->getParameter());
466 std::vector<Real> Bz(FEM_->nu(),0.0);
467 FEM_->apply_jacobian_2(Bz,*zp,this->getParameter());
468 plus(F,Bz);
469 // Solve solve state
470 FEM_->linear_solve(*up,dl,d,du,F);
471 num_solves_++;
472 // Compute residual
473 value(c,u,z,tol);
474 }
475
477 const ROL::Vector<Real> &z, Real &tol) {
478 ROL::Ptr<std::vector<Real> > jvp = dynamic_cast<ROL::StdVector<Real>&>(jv).getVector();
479 ROL::Ptr<const std::vector<Real> > vp = dynamic_cast<const ROL::StdVector<Real>&>(v).getVector();
480 FEM_->apply_jacobian_1(*jvp, *vp, this->getParameter());
481 }
482
484 const ROL::Vector<Real> &z, Real &tol) {
485 ROL::Ptr<std::vector<Real> > jvp = dynamic_cast<ROL::StdVector<Real>&>(jv).getVector();
486 ROL::Ptr<const std::vector<Real> > vp = dynamic_cast<const ROL::StdVector<Real>&>(v).getVector();
487 FEM_->apply_jacobian_2(*jvp, *vp, this->getParameter());
488 scale(*jvp,-1.0);
489 }
490
492 const ROL::Vector<Real> &z, Real &tol) {
493 ROL::Ptr<std::vector<Real> > jvp = dynamic_cast<ROL::StdVector<Real>&>(jv).getVector();
494 ROL::Ptr<const std::vector<Real> > vp = dynamic_cast<const ROL::StdVector<Real>&>(v).getVector();
495 // Get PDE Jacobian
496 std::vector<Real> dl(FEM_->nu()-1,0.0);
497 std::vector<Real> d(FEM_->nu(),0.0);
498 std::vector<Real> du(FEM_->nu()-1,0.0);
499 FEM_->build_jacobian_1(dl,d,du,this->getParameter());
500 // Solve solve state
501 FEM_->linear_solve(*jvp,dl,d,du,*vp);
502 num_solves_++;
503 }
504
506 const ROL::Vector<Real> &z, Real &tol) {
507 applyJacobian_1(jv,v,u,z,tol);
508 }
509
511 const ROL::Vector<Real> &z, Real &tol) {
512 ROL::Ptr<std::vector<Real> > jvp = dynamic_cast<ROL::StdVector<Real>&>(jv).getVector();
513 ROL::Ptr<const std::vector<Real> > vp = dynamic_cast<const ROL::StdVector<Real>&>(v).getVector();
514 FEM_->apply_adjoint_jacobian_2(*jvp, *vp, this->getParameter());
515 scale(*jvp,-1.0);
516 }
517
519 const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
520 applyInverseJacobian_1(jv,v,u,z,tol);
521 }
522
524 const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
525 ahwv.zero();
526 }
527
529 const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
530 ahwv.zero();
531 }
532
534 const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
535 ahwv.zero();
536 }
537
539 const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
540 ahwv.zero();
541 }
542};
543
544// Class BurgersObjective contains the necessary information
545// to compute the value and gradient of the objective function.
546template<class Real>
548private:
549 const ROL::Ptr<FEM<Real> > FEM_;
550 const Real alpha_;
551
552/***************************************************************/
553/********** BEGIN PRIVATE MEMBER FUNCTION DECLARATION **********/
554/***************************************************************/
555 // Evaluates the discretized L2 inner product of two vectors.
556 Real dot(const std::vector<Real> &x, const std::vector<Real> &y) {
557 Real ip = 0.0;
558 int size = static_cast<int>(x.size());
559 for (int i=0; i<size; i++) {
560 ip += x[i]*y[i];
561 }
562 return ip;
563 }
564
565 // Returns u <- alpha*u where u is a vector and alpha is a scalar
566 void scale(std::vector<Real> &u, const Real alpha=0.0) {
567 int size = static_cast<int>(u.size());
568 for (int i=0; i<size; i++) {
569 u[i] *= alpha;
570 }
571 }
572/*************************************************************/
573/********** END PRIVATE MEMBER FUNCTION DECLARATION **********/
574/*************************************************************/
575
576public:
577 DiffusionObjective(const ROL::Ptr<FEM<Real> > fem, const Real alpha = 1.e-4)
578 : FEM_(fem), alpha_(alpha) {}
579
580 Real value( const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
581 ROL::Ptr<const std::vector<Real> > up = dynamic_cast<const ROL::StdVector<Real>&>(u).getVector();
582 ROL::Ptr<const std::vector<Real> > zp = dynamic_cast<const ROL::StdVector<Real>&>(z).getVector();
583 int nz = FEM_->nz(), nu = FEM_->nu();
584 // COMPUTE CONTROL PENALTY
585 std::vector<Real> Mz(nz);
586 FEM_->apply_mass_control(Mz,*zp);
587 Real zval = alpha_*0.5*dot(Mz,*zp);
588 // COMPUTE STATE TRACKING TERM
589 std::vector<Real> diff(nu);
590 for (int i=0; i<nu; i++) {
591 diff[i] = ((*up)[i]-FEM_->evaluate_target(i+1,this->getParameter()));
592 }
593 std::vector<Real> Mu(nu);
594 FEM_->apply_mass_state(Mu,diff,this->getParameter());
595 Real uval = 0.5*dot(Mu,diff);
596 return uval+zval;
597 }
598
599 void gradient_1( ROL::Vector<Real> &g, const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
600 ROL::Ptr<std::vector<Real> > gp = dynamic_cast<ROL::StdVector<Real>&>(g).getVector();
601 ROL::Ptr<const std::vector<Real> > up = dynamic_cast<const ROL::StdVector<Real>&>(u).getVector();
602 int nu = FEM_->nu();
603 // COMPUTE GRADIENT
604 std::vector<Real> diff(nu);
605 for (int i=0; i<nu; i++) {
606 diff[i] = ((*up)[i]-FEM_->evaluate_target(i+1,this->getParameter()));
607 }
608 FEM_->apply_mass_state(*gp,diff,this->getParameter());
609 }
610
611 void gradient_2( ROL::Vector<Real> &g, const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
612 ROL::Ptr<std::vector<Real> > gp = dynamic_cast<ROL::StdVector<Real>&>(g).getVector();
613 ROL::Ptr<const std::vector<Real> > zp = dynamic_cast<const ROL::StdVector<Real>&>(z).getVector();
614 // COMPUTE GRADIENT
615 FEM_->apply_mass_control(*gp,*zp);
616 scale(*gp,alpha_);
617 }
618
620 const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
621 ROL::Ptr<std::vector<Real> > hvp = dynamic_cast<ROL::StdVector<Real>&>(hv).getVector();
622 ROL::Ptr<const std::vector<Real> > vp = dynamic_cast<const ROL::StdVector<Real>&>(v).getVector();
623 // COMPUTE HESSVEC
624 FEM_->apply_mass_state(*hvp,*vp,this->getParameter());
625 }
626
628 const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
629 hv.zero();
630 }
631
633 const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
634 hv.zero();
635 }
636
638 const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol ) {
639 ROL::Ptr<std::vector<Real> > hvp = dynamic_cast<ROL::StdVector<Real>&>(hv).getVector();
640 ROL::Ptr<const std::vector<Real> > vp = dynamic_cast<const ROL::StdVector<Real>&>(v).getVector();
641 // COMPUTE HESSVEC
642 FEM_->apply_mass_control(*hvp,*vp);
643 scale(*hvp,alpha_);
644 }
645
646};
Contains definitions of custom data types in ROL.
void applyAdjointHessian_11(ROL::Vector< Real > &ahwv, const ROL::Vector< Real > &w, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Apply the simulation-space derivative of the adjoint of the constraint simulation-space Jacobian at ...
void applyAdjointHessian_12(ROL::Vector< Real > &ahwv, const ROL::Vector< Real > &w, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Apply the optimization-space derivative of the adjoint of the constraint simulation-space Jacobian at...
const ROL::Ptr< FEM< Real > > FEM_
void applyJacobian_2(ROL::Vector< Real > &jv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Apply the partial constraint Jacobian at , , to the vector .
DiffusionConstraint(const ROL::Ptr< FEM< Real > > &FEM)
void applyJacobian_1(ROL::Vector< Real > &jv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Apply the partial constraint Jacobian at , , to the vector .
void scale(std::vector< Real > &u, const Real alpha=0.0)
void applyInverseAdjointJacobian_1(ROL::Vector< Real > &jv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Apply the inverse of the adjoint of the partial constraint Jacobian at , , to the vector .
void applyInverseJacobian_1(ROL::Vector< Real > &jv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Apply the inverse partial constraint Jacobian at , , to the vector .
void value(ROL::Vector< Real > &c, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Evaluate the constraint operator at .
void applyAdjointJacobian_1(ROL::Vector< Real > &jv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Apply the adjoint of the partial constraint Jacobian at , , to the vector . This is the primary inter...
void applyAdjointHessian_21(ROL::Vector< Real > &ahwv, const ROL::Vector< Real > &w, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Apply the simulation-space derivative of the adjoint of the constraint optimization-space Jacobian at...
void plus(std::vector< Real > &u, const std::vector< Real > &s, const Real alpha=1.0)
void applyAdjointHessian_22(ROL::Vector< Real > &ahwv, const ROL::Vector< Real > &w, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Apply the optimization-space derivative of the adjoint of the constraint optimization-space Jacobian ...
void solve(ROL::Vector< Real > &c, ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Given , solve for .
void applyAdjointJacobian_2(ROL::Vector< Real > &jv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Apply the adjoint of the partial constraint Jacobian at , , to vector . This is the primary interface...
void gradient_2(ROL::Vector< Real > &g, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Compute gradient with respect to second component.
const ROL::Ptr< FEM< Real > > FEM_
void hessVec_11(ROL::Vector< Real > &hv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Apply Hessian approximation to vector.
void gradient_1(ROL::Vector< Real > &g, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Compute gradient with respect to first component.
void hessVec_21(ROL::Vector< Real > &hv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
void hessVec_22(ROL::Vector< Real > &hv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
DiffusionObjective(const ROL::Ptr< FEM< Real > > fem, const Real alpha=1.e-4)
Real value(const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Compute value.
void hessVec_12(ROL::Vector< Real > &hv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Real dot(const std::vector< Real > &x, const std::vector< Real > &y)
void scale(std::vector< Real > &u, const Real alpha=0.0)
void build_mesh(std::vector< Real > &x, const std::vector< Real > &param)
std::vector< Real > pts_
void build_mesh(std::vector< Real > &x)
void apply_adjoint_jacobian_2(std::vector< Real > &jv, const std::vector< Real > &v, const std::vector< Real > &param)
void apply_jacobian_2(std::vector< Real > &jv, const std::vector< Real > &v, const std::vector< Real > &param)
FEM(int nx=128, Real kl=0.1, Real kr=10.0)
void apply_jacobian_1(std::vector< Real > &jv, const std::vector< Real > &v, const std::vector< Real > &param)
void apply_inverse_mass_control(std::vector< Real > &Mv, const std::vector< Real > &v)
Real evaluate_target(int i, const std::vector< Real > &param)
void apply_mass_state(std::vector< Real > &Mv, const std::vector< Real > &v, const std::vector< Real > &param)
void build_jacobian_1(std::vector< Real > &dl, std::vector< Real > &d, std::vector< Real > &du, const std::vector< Real > &param)
std::vector< Real > wts_
void build_force(std::vector< Real > &F, const std::vector< Real > &param)
void linear_solve(std::vector< Real > &u, std::vector< Real > &dl, std::vector< Real > &d, std::vector< Real > &du, const std::vector< Real > &r, const bool transpose=false)
void apply_mass_control(std::vector< Real > &Mv, const std::vector< Real > &v)
Defines the constraint operator interface for simulation-based optimization.
const std::vector< Real > getParameter(void) const
Provides the interface to evaluate simulation-based objective functions.
const std::vector< Real > getParameter(void) const
Provides the ROL::Vector interface for scalar values, to be used, for example, with scalar constraint...
Defines the linear algebra or vector space interface.
virtual void zero()
Set to zero vector.