58 #include "Teuchos_oblackholestream.hpp" 59 #include "Teuchos_GlobalMPISession.hpp" 60 #include "Teuchos_XMLParameterListHelpers.hpp" 61 #include "Teuchos_LAPACK.hpp" 85 return std::sqrt(
dot(r,r));
88 Real
dot(
const std::vector<Real> &x,
const std::vector<Real> &y) {
90 Real c = ((x.size()==
nx_) ? 4.0 : 2.0);
91 for (
unsigned i = 0; i < x.size(); i++) {
93 ip += dx_/6.0*(c*x[i] + x[i+1])*y[i];
95 else if ( i == x.size()-1 ) {
96 ip += dx_/6.0*(x[i-1] + c*x[i])*y[i];
99 ip += dx_/6.0*(x[i-1] + 4.0*x[i] + x[i+1])*y[i];
105 void update(std::vector<Real> &u,
const std::vector<Real> &s,
const Real alpha=1.0) {
106 for (
unsigned i = 0; i < u.size(); i++) {
111 void scale(std::vector<Real> &u,
const Real alpha=0.0) {
112 for (
unsigned i = 0; i < u.size(); i++) {
117 void compute_residual(std::vector<Real> &r,
const std::vector<Real> &uold,
const std::vector<Real> &zold,
118 const std::vector<Real> &unew,
const std::vector<Real> &znew) {
121 for (
unsigned n = 0; n <
nx_; n++) {
123 r[n] += (4.0*dx_/6.0 + 0.5*dt_*2.0*nu_/
dx_)*unew[n];
124 r[n] += (-4.0*dx_/6.0 + 0.5*dt_*2.0*nu_/
dx_)*uold[n];
126 r[n] += (dx_/6.0 - 0.5*dt_*nu_/
dx_)*unew[n-1];
127 r[n] -= (dx_/6.0 + 0.5*dt_*nu_/
dx_)*uold[n-1];
130 r[n] += (dx_/6.0 - 0.5*dt_*nu_/
dx_)*unew[n+1];
131 r[n] -= (dx_/6.0 + 0.5*dt_*nu_/
dx_)*uold[n+1];
135 r[n] -= 0.5*dt_*unew[n-1]*(unew[n-1]+unew[n])/6.0;
136 r[n] -= 0.5*dt_*uold[n-1]*(uold[n-1]+uold[n])/6.0;
139 r[n] += 0.5*dt_*unew[n+1]*(unew[n]+unew[n+1])/6.0;
140 r[n] += 0.5*dt_*uold[n+1]*(uold[n]+uold[n+1])/6.0;
143 r[n] -= 0.5*dt_*dx_/6.0*(znew[n]+4.0*znew[n+1]+znew[n+2]);
144 r[n] -= 0.5*dt_*dx_/6.0*(zold[n]+4.0*zold[n+1]+zold[n+2]);
149 r[0] -= dt_*(0.5*u0_*(unew[ 0] + uold[ 0])/6.0 + u0_*u0_/6.0 + nu_*u0_/dx_);
150 r[nx_-1] += dt_*(0.5*u1_*(unew[nx_-1] + uold[nx_-1])/6.0 + u1_*u1_/6.0 - nu_*u1_/dx_);
154 const std::vector<Real> &u) {
157 d.resize(nx_,4.0*dx_/6.0 + 0.5*dt_*nu_*2.0/dx_);
159 dl.resize(nx_-1,dx_/6.0 - 0.5*dt_*nu_/dx_);
161 du.resize(nx_-1,dx_/6.0 - 0.5*dt_*nu_/dx_);
163 for (
unsigned n = 0; n <
nx_; n++) {
165 dl[n] += 0.5*dt_*(-2.0*u[n]-u[n+1])/6.0;
166 d[n] += 0.5*dt_*u[n+1]/6.0;
169 d[n] -= 0.5*dt_*u[n-1]/6.0;
170 du[n-1] += 0.5*dt_*(u[n-1]+2.0*u[n])/6.0;
174 d[0] -= 0.5*dt_*u0_/6.0;
175 d[nx_-1] += 0.5*dt_*u1_/6.0;
179 bool adjoint =
false) {
183 for (
unsigned n = 0; n <
nx_; n++) {
184 jv[n] = (4.0*dx_/6.0 + 0.5*dt_*nu_/dx_*2.0)*v[n];
186 jv[n] += dx_/6.0*v[n-1]-0.5*dt_*nu_/dx_*v[n-1];
188 jv[n] -= 0.5*dt_*(u[n-1]/6.0*v[n]-(u[n-1]+2.0*u[n])/6.0*v[n-1]);
191 jv[n] -= 0.5*dt_*(u[n-1]/6.0*v[n]+(u[n]+2.0*u[n-1])/6.0*v[n-1]);
195 jv[n] += dx_/6.0*v[n+1]-0.5*dt_*nu_/dx_*v[n+1];
197 jv[n] += 0.5*dt_*(u[n+1]/6.0*v[n]-(u[n+1]+2.0*u[n])/6.0*v[n+1]);
200 jv[n] += 0.5*dt_*(u[n+1]/6.0*v[n]+(u[n]+2.0*u[n+1])/6.0*v[n+1]);
204 jv[0] -= 0.5*dt_*u0_/6.0*v[0];
205 jv[nx_-1] += 0.5*dt_*u1_/6.0*v[nx_-1];
209 bool adjoint =
false) {
213 for (
unsigned n = 0; n <
nx_; n++) {
214 jv[n] = (-4.0*dx_/6.0 + 0.5*dt_*nu_/dx_*2.0)*v[n];
216 jv[n] += -dx_/6.0*v[n-1]-0.5*dt_*nu_/dx_*v[n-1];
218 jv[n] -= 0.5*dt_*(u[n-1]/6.0*v[n]-(u[n-1]+2.0*u[n])/6.0*v[n-1]);
221 jv[n] -= 0.5*dt_*(u[n-1]/6.0*v[n]+(u[n]+2.0*u[n-1])/6.0*v[n-1]);
225 jv[n] += -dx_/6.0*v[n+1]-0.5*dt_*nu_/dx_*v[n+1];
227 jv[n] += 0.5*dt_*(u[n+1]/6.0*v[n]-(u[n+1]+2.0*u[n])/6.0*v[n+1]);
230 jv[n] += 0.5*dt_*(u[n+1]/6.0*v[n]+(u[n]+2.0*u[n+1])/6.0*v[n+1]);
234 jv[0] -= 0.5*dt_*u0_/6.0*v[0];
235 jv[nx_-1] += 0.5*dt_*u1_/6.0*v[nx_-1];
238 void apply_pde_jacobian(std::vector<Real> &jv,
const std::vector<Real> &vold,
const std::vector<Real> &uold,
239 const std::vector<Real> &vnew,
const std::vector<Real> unew,
bool adjoint =
false) {
243 for (
unsigned n = 0; n <
nx_; n++) {
244 jv[n] += (4.0*dx_/6.0+0.5*dt_*nu_/dx_*2.0)*vnew[n];
245 jv[n] -= (4.0*dx_/6.0-0.5*dt_*nu_/dx_*2.0)*vold[n];
247 jv[n] += dx_/6.0*vnew[n-1]-0.5*dt_*nu_/dx_*vnew[n-1];
248 jv[n] -= dx_/6.0*vold[n-1]+0.5*dt_*nu_/dx_*vold[n-1];
250 jv[n] -= 0.5*dt_*(unew[n-1]/6.0*vnew[n]-(unew[n-1]+2.0*unew[n])/6.0*vnew[n-1]);
251 jv[n] -= 0.5*dt_*(uold[n-1]/6.0*vold[n]-(uold[n-1]+2.0*uold[n])/6.0*vold[n-1]);
254 jv[n] -= 0.5*dt_*(unew[n-1]/6.0*vnew[n]+(unew[n]+2.0*unew[n-1])/6.0*vnew[n-1]);
255 jv[n] -= 0.5*dt_*(uold[n-1]/6.0*vold[n]+(uold[n]+2.0*uold[n-1])/6.0*vold[n-1]);
259 jv[n] += dx_/6.0*vnew[n+1]-0.5*dt_*nu_/dx_*vnew[n+1];
260 jv[n] -= dx_/6.0*vold[n+1]+0.5*dt_*nu_/dx_*vold[n+1];
262 jv[n] += 0.5*dt_*(unew[n+1]/6.0*vnew[n]-(unew[n+1]+2.0*unew[n])/6.0*vnew[n+1]);
263 jv[n] += 0.5*dt_*(uold[n+1]/6.0*vold[n]-(uold[n+1]+2.0*uold[n])/6.0*vold[n+1]);
266 jv[n] += 0.5*dt_*(unew[n+1]/6.0*vnew[n]+(unew[n]+2.0*unew[n+1])/6.0*vnew[n+1]);
267 jv[n] += 0.5*dt_*(uold[n+1]/6.0*vold[n]+(uold[n]+2.0*uold[n+1])/6.0*vold[n+1]);
271 jv[0] -= 0.5*dt_*u0_/6.0*vnew[0];
272 jv[0] -= 0.5*dt_*u0_/6.0*vold[0];
273 jv[nx_-1] += 0.5*dt_*u1_/6.0*vnew[nx_-1];
274 jv[nx_-1] += 0.5*dt_*u1_/6.0*vold[nx_-1];
277 void apply_pde_hessian(std::vector<Real> &hv,
const std::vector<Real> &wold,
const std::vector<Real> &vold,
278 const std::vector<Real> &wnew,
const std::vector<Real> &vnew) {
281 for (
unsigned n = 0; n <
nx_; n++) {
283 hv[n] += 0.5*dt_*((wnew[n-1]*(vnew[n-1]+2.0*vnew[n]) - wnew[n]*vnew[n-1])/6.0);
284 hv[n] += 0.5*dt_*((wold[n-1]*(vold[n-1]+2.0*vold[n]) - wold[n]*vold[n-1])/6.0);
287 hv[n] += 0.5*dt_*((wnew[n]*vnew[n+1] - wnew[n+1]*(2.0*vnew[n]+vnew[n+1]))/6.0);
288 hv[n] += 0.5*dt_*((wold[n]*vold[n+1] - wold[n+1]*(2.0*vold[n]+vold[n+1]))/6.0);
295 unsigned dim = ((adjoint ==
true) ? nx_+2 : nx_);
297 for (
unsigned n = 0; n < dim; n++) {
300 jv[n] = -0.5*dt_*dx_/6.0*v[n];
303 jv[n] = -0.5*dt_*dx_/6.0*(4.0*v[n-1]+v[n]);
305 else if ( n == nx_ ) {
306 jv[n] = -0.5*dt_*dx_/6.0*(4.0*v[n-1]+v[n-2]);
308 else if ( n == nx_+1 ) {
309 jv[n] = -0.5*dt_*dx_/6.0*v[n-2];
312 jv[n] = -0.5*dt_*dx_/6.0*(v[n-2]+4.0*v[n-1]+v[n]);
316 jv[n] -= 0.5*dt_*dx_/6.0*(v[n]+4.0*v[n+1]+v[n+2]);
321 void run_newton(std::vector<Real> &u,
const std::vector<Real> &znew,
322 const std::vector<Real> &uold,
const std::vector<Real> &zold) {
326 std::vector<Real> r(nx_,0.0);
331 unsigned maxit = 500;
333 std::vector<Real> d(nx_,0.0);
334 std::vector<Real> dl(nx_-1,0.0);
335 std::vector<Real> du(nx_-1,0.0);
337 Real alpha = 1.0, tmp = 0.0;
338 std::vector<Real> s(nx_,0.0);
339 std::vector<Real> utmp(nx_,0.0);
340 for (
unsigned i = 0; i < maxit; i++) {
349 utmp.assign(u.begin(),u.end());
353 while ( rnorm > (1.0-1.e-4*alpha)*tmp && alpha > std::sqrt(
ROL::ROL_EPSILON) ) {
355 utmp.assign(u.begin(),u.end());
361 u.assign(utmp.begin(),utmp.end());
362 if ( rnorm < rtol ) {
369 const std::vector<Real> &dl,
const std::vector<Real> &d,
const std::vector<Real> &du,
370 const std::vector<Real> &r,
const bool transpose =
false) {
371 bool useLAPACK =
false;
373 u.assign(r.begin(),r.end());
375 std::vector<Real> Dl(dl);
376 std::vector<Real> Du(du);
377 std::vector<Real> D(d);
379 Teuchos::LAPACK<int,Real> lp;
380 std::vector<Real> Du2(nx_-2,0.0);
381 std::vector<int> ipiv(nx_,0);
385 lp.GTTRF(nx_,&Dl[0],&D[0],&Du[0],&Du2[0],&ipiv[0],&info);
386 char trans = ((transpose ==
true) ?
'T' :
'N');
387 lp.GTTRS(trans,nx_,nhrs,&Dl[0],&D[0],&Du[0],&Du2[0],&ipiv[0],&u[0],ldb,&info);
392 unsigned maxit = 100;
393 Real rtol = std::min(1.e-12,1.e-4*std::sqrt(
dot(r,r)));
395 Real rnorm = 10.0*rtol;
396 for (
unsigned i = 0; i < maxit; i++) {
397 for (
unsigned n = 0; n <
nx_; n++) {
400 u[n] -= ((transpose ==
false) ? du[n] : dl[n])*u[n+1]/d[n];
403 u[n] -= ((transpose ==
false) ? dl[n-1] : du[n-1])*u[n-1]/d[n];
408 for (
unsigned n = 0; n <
nx_; n++) {
409 resid = r[n] - d[n]*u[n];
411 resid -= ((transpose ==
false) ? du[n] : dl[n])*u[n+1];
414 resid -= ((transpose ==
false) ? dl[n-1] : du[n-1])*u[n-1];
416 rnorm += resid*resid;
418 rnorm = std::sqrt(rnorm);
419 if ( rnorm < rtol ) {
430 Real nu = 1.e-2, Real u0 = 0.0, Real u1 = 0.0, Real f = 0.0)
431 : nx_((unsigned)nx), nt_((unsigned)nt), T_(T), nu_(nu), u0_(u0), u1_(u1), f_(f) {
432 dx_ = 1.0/((Real)nx+1.0);
435 u_init_.resize(nx_,0.0);
437 for (
unsigned n = 0; n <
nx_; n++) {
439 u_init_[n] = ((x <= 0.5) ? 1.0 : 0.0);
444 Teuchos::RCP<std::vector<Real> > cp =
445 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<
ROL::StdVector<Real> >(c)).getVector());
446 Teuchos::RCP<const std::vector<Real> > up =
448 Teuchos::RCP<const std::vector<Real> > zp =
451 std::vector<Real> C(nx_,0.0);
452 std::vector<Real> uold(u_init_);
453 std::vector<Real> unew(u_init_);
454 std::vector<Real> zold(nx_+2,0.0);
455 std::vector<Real> znew(nx_+2,0.0);
457 for (
unsigned n = 0; n < nx_+2; n++) {
461 for (
unsigned t = 0; t <
nt_; t++) {
463 for (
unsigned n = 0; n <
nx_; n++) {
464 unew[n] = (*up)[t*nx_+n];
467 for (
unsigned n = 0; n < nx_+2; n++) {
468 znew[n] = (*zp)[(t+1)*(nx_+2)+n];
473 for (
unsigned n = 0; n <
nx_; n++) {
474 (*cp)[t*nx_+n] = C[n];
477 uold.assign(unew.begin(),unew.end());
478 zold.assign(znew.begin(),znew.end());
483 Teuchos::RCP<std::vector<Real> > up =
484 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<
ROL::StdVector<Real> >(u)).getVector());
485 up->assign(up->size(),z.
norm()/up->size());
486 Teuchos::RCP<const std::vector<Real> > zp =
489 std::vector<Real> uold(u_init_);
490 std::vector<Real> unew(u_init_);
491 std::vector<Real> zold(nx_+2,0.0);
492 std::vector<Real> znew(nx_+2,0.0);
494 for (
unsigned n = 0; n < nx_+2; n++) {
498 for (
unsigned t = 0; t <
nt_; t++) {
500 for (
unsigned n = 0; n < nx_+2; n++) {
501 znew[n] = (*zp)[(t+1)*(nx_+2)+n];
506 for (
unsigned n = 0; n <
nx_; n++) {
507 (*up)[t*nx_+n] = unew[n];
510 uold.assign(unew.begin(),unew.end());
511 zold.assign(znew.begin(),znew.end());
518 Teuchos::RCP<std::vector<Real> > jvp =
519 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<
ROL::StdVector<Real> >(jv)).getVector());
520 Teuchos::RCP<const std::vector<Real> > vp =
522 Teuchos::RCP<const std::vector<Real> > up =
524 std::vector<Real> J(nx_,0.0);
525 std::vector<Real> vold(nx_,0.0);
526 std::vector<Real> vnew(nx_,0.0);
527 std::vector<Real> uold(nx_,0.0);
528 std::vector<Real> unew(nx_,0.0);
529 for (
unsigned t = 0; t <
nt_; t++) {
530 for (
unsigned n = 0; n <
nx_; n++) {
531 unew[n] = (*up)[t*nx_+n];
532 vnew[n] = (*vp)[t*nx_+n];
535 for (
unsigned n = 0; n <
nx_; n++) {
536 (*jvp)[t*nx_+n] = J[n];
538 vold.assign(vnew.begin(),vnew.end());
539 uold.assign(unew.begin(),unew.end());
545 Teuchos::RCP<std::vector<Real> > jvp =
546 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<
ROL::StdVector<Real> >(jv)).getVector());
547 Teuchos::RCP<const std::vector<Real> > vp =
549 Teuchos::RCP<const std::vector<Real> > zp =
551 std::vector<Real> vnew(nx_+2,0.0);
552 std::vector<Real> vold(nx_+2,0.0);
553 std::vector<Real> jnew(nx_,0.0);
554 std::vector<Real> jold(nx_,0.0);
555 for (
unsigned n = 0; n < nx_+2; n++) {
559 for (
unsigned t = 0; t <
nt_; t++) {
560 for (
unsigned n = 0; n < nx_+2; n++) {
561 vnew[n] = (*vp)[(t+1)*(nx_+2)+n];
564 for (
unsigned n = 0; n <
nx_; n++) {
566 (*jvp)[t*nx_+n] = jnew[n] + jold[n];
568 jold.assign(jnew.begin(),jnew.end());
574 Teuchos::RCP<std::vector<Real> > ijvp =
575 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<
ROL::StdVector<Real> >(ijv)).getVector());
576 Teuchos::RCP<const std::vector<Real> > vp =
578 Teuchos::RCP<const std::vector<Real> > up =
580 std::vector<Real> J(nx_,0.0);
581 std::vector<Real> r(nx_,0.0);
582 std::vector<Real> s(nx_,0.0);
583 std::vector<Real> uold(nx_,0.0);
584 std::vector<Real> unew(nx_,0.0);
585 std::vector<Real> d(nx_,0.0);
586 std::vector<Real> dl(nx_-1,0.0);
587 std::vector<Real> du(nx_-1,0.0);
588 for (
unsigned t = 0; t <
nt_; t++) {
590 for (
unsigned n = 0; n <
nx_; n++) {
591 r[n] = (*vp)[t*nx_+n];
592 unew[n] = (*up)[t*nx_+n];
601 for (
unsigned n = 0; n <
nx_; n++) {
602 (*ijvp)[t*nx_+n] = s[n];
604 uold.assign(unew.begin(),unew.end());
610 Teuchos::RCP<std::vector<Real> > jvp =
611 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<
ROL::StdVector<Real> >(ajv)).getVector());
612 Teuchos::RCP<const std::vector<Real> > vp =
614 Teuchos::RCP<const std::vector<Real> > up =
616 std::vector<Real> J(nx_,0.0);
617 std::vector<Real> vold(nx_,0.0);
618 std::vector<Real> vnew(nx_,0.0);
619 std::vector<Real> unew(nx_,0.0);
620 for (
unsigned t = nt_; t > 0; t--) {
621 for (
unsigned n = 0; n <
nx_; n++) {
622 unew[n] = (*up)[(t-1)*nx_+n];
623 vnew[n] = (*vp)[(t-1)*nx_+n];
626 for (
unsigned n = 0; n <
nx_; n++) {
627 (*jvp)[(t-1)*nx_+n] = J[n];
629 vold.assign(vnew.begin(),vnew.end());
635 Teuchos::RCP<std::vector<Real> > jvp =
636 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<
ROL::StdVector<Real> >(jv)).getVector());
637 Teuchos::RCP<const std::vector<Real> > vp =
639 Teuchos::RCP<const std::vector<Real> > zp =
641 std::vector<Real> vnew(nx_,0.0);
642 std::vector<Real> vold(nx_,0.0);
643 std::vector<Real> jnew(nx_+2,0.0);
644 std::vector<Real> jold(nx_+2,0.0);
645 for (
unsigned t = nt_+1; t > 0; t--) {
646 for (
unsigned n = 0; n <
nx_; n++) {
648 vnew[n] = (*vp)[(t-2)*nx_+n];
655 for (
unsigned n = 0; n < nx_+2; n++) {
657 (*jvp)[(t-1)*(nx_+2)+n] = jnew[n] + jold[n];
659 jold.assign(jnew.begin(),jnew.end());
665 Teuchos::RCP<std::vector<Real> > ijvp =
666 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<
ROL::StdVector<Real> >(ijv)).getVector());
667 Teuchos::RCP<const std::vector<Real> > vp =
669 Teuchos::RCP<const std::vector<Real> > up =
671 std::vector<Real> J(nx_,0.0);
672 std::vector<Real> r(nx_,0.0);
673 std::vector<Real> s(nx_,0.0);
674 std::vector<Real> unew(nx_,0.0);
675 std::vector<Real> d(nx_,0.0);
676 std::vector<Real> dl(nx_-1,0.0);
677 std::vector<Real> du(nx_-1,0.0);
678 for (
unsigned t = nt_; t > 0; t--) {
680 for (
unsigned n = 0; n <
nx_; n++) {
681 r[n] = (*vp)[(t-1)*nx_+n];
682 unew[n] = (*up)[(t-1)*nx_+n];
691 for (
unsigned n = 0; n <
nx_; n++) {
692 (*ijvp)[(t-1)*nx_+n] = s[n];
699 Teuchos::RCP<std::vector<Real> > hwvp =
700 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<
ROL::StdVector<Real> >(hwv)).getVector());
701 Teuchos::RCP<const std::vector<Real> > wp =
703 Teuchos::RCP<const std::vector<Real> > vp =
705 std::vector<Real> snew(nx_,0.0);
706 std::vector<Real> wnew(nx_,0.0);
707 std::vector<Real> wold(nx_,0.0);
708 std::vector<Real> vnew(nx_,0.0);
709 for (
unsigned t = nt_; t > 0; t--) {
710 for (
unsigned n = 0; n <
nx_; n++) {
711 vnew[n] = (*vp)[(t-1)*nx_+n];
712 wnew[n] = (*wp)[(t-1)*nx_+n];
715 for (
unsigned n = 0; n <
nx_; n++) {
716 (*hwvp)[(t-1)*nx_+n] = snew[n];
718 wold.assign(wnew.begin(),wnew.end());
754 case 1: val = ((x<=0.5) ? 1.0 : 0.0);
break;
755 case 2: val = 1.0;
break;
756 case 3: val = std::abs(std::sin(8.0*M_PI*x));
break;
757 case 4: val = std::exp(-0.5*(x-0.5)*(x-0.5));
break;
762 Real
dot(
const std::vector<Real> &x,
const std::vector<Real> &y) {
764 Real c = ((x.size()==
nx_) ? 4.0 : 2.0);
765 for (
unsigned i=0; i<x.size(); i++) {
767 ip += dx_/6.0*(c*x[i] + x[i+1])*y[i];
769 else if ( i == x.size()-1 ) {
770 ip += dx_/6.0*(x[i-1] + c*x[i])*y[i];
773 ip += dx_/6.0*(x[i-1] + 4.0*x[i] + x[i+1])*y[i];
779 void apply_mass(std::vector<Real> &Mu,
const std::vector<Real> &u ) {
780 Mu.resize(u.size(),0.0);
781 Real c = ((u.size()==
nx_) ? 4.0 : 2.0);
782 for (
unsigned i=0; i<u.size(); i++) {
784 Mu[i] = dx_/6.0*(c*u[i] + u[i+1]);
786 else if ( i == u.size()-1 ) {
787 Mu[i] = dx_/6.0*(u[i-1] + c*u[i]);
790 Mu[i] = dx_/6.0*(u[i-1] + 4.0*u[i] + u[i+1]);
801 : alpha_(alpha), nx_((unsigned)nx), nt_((unsigned)nt), T_(T) {
802 dx_ = 1.0/((Real)nx+1.0);
807 Teuchos::RCP<const std::vector<Real> > up =
809 Teuchos::RCP<const std::vector<Real> > zp =
812 std::vector<Real> U(nx_,0.0);
813 std::vector<Real> Z(nx_+2,0.0);
814 for (
unsigned n = 0; n < nx_+2; n++) {
818 Real val = 0.5*ss*alpha_*
dot(Z,Z);
819 for (
unsigned t = 0; t <
nt_; t++) {
820 ss = ((t < nt_-1) ? dt_ : 0.5*dt_);
821 for (
unsigned n = 0; n <
nx_; n++) {
822 U[n] = (*up)[t*nx_+n]-evaluate_target((Real)(n+1)*dx_);
824 val += 0.5*ss*
dot(U,U);
825 for (
unsigned n = 0; n < nx_+2; n++) {
826 Z[n] = (*zp)[(t+1)*(nx_+2)+n];
828 val += 0.5*ss*alpha_*
dot(Z,Z);
834 Teuchos::RCP<std::vector<Real> > gp = Teuchos::rcp_const_cast<std::vector<Real> >(
836 Teuchos::RCP<const std::vector<Real> > up =
839 std::vector<Real> U(nx_,0.0);
840 std::vector<Real> M(nx_,0.0);
842 for (
unsigned t = 0; t <
nt_; t++) {
843 ss = ((t < nt_-1) ? dt_ : 0.5*dt_);
844 for (
unsigned n = 0; n <
nx_; n++) {
845 U[n] = (*up)[t*nx_+n]-evaluate_target((Real)(n+1)*dx_);
848 for (
unsigned n = 0; n <
nx_; n++) {
849 (*gp)[t*nx_+n] = ss*M[n];
855 Teuchos::RCP<std::vector<Real> > gp = Teuchos::rcp_const_cast<std::vector<Real> >(
857 Teuchos::RCP<const std::vector<Real> > zp =
860 std::vector<Real> Z(nx_+2,0.0);
861 std::vector<Real> M(nx_+2,0.0);
863 for (
unsigned t = 0; t < nt_+1; t++) {
864 ss = ((t < nt_ && t > 0) ? dt_ : 0.5*dt_);
865 for (
unsigned n = 0; n < nx_+2; n++) {
866 Z[n] = (*zp)[t*(nx_+2)+n];
869 for (
unsigned n = 0; n < nx_+2; n++) {
870 (*gp)[t*(nx_+2)+n] = ss*alpha_*M[n];
877 Teuchos::RCP<std::vector<Real> > hvp = Teuchos::rcp_const_cast<std::vector<Real> >(
879 Teuchos::RCP<const std::vector<Real> > vp =
882 std::vector<Real> V(nx_,0.0);
883 std::vector<Real> M(nx_,0.0);
885 for (
unsigned t = 0; t <
nt_; t++) {
886 ss = ((t < nt_-1) ? dt_ : 0.5*dt_);
887 for (
unsigned n = 0; n <
nx_; n++) {
888 V[n] = (*vp)[t*nx_+n];
891 for (
unsigned n = 0; n <
nx_; n++) {
892 (*hvp)[t*nx_+n] = ss*M[n];
909 Teuchos::RCP<std::vector<Real> > hvp = Teuchos::rcp_const_cast<std::vector<Real> >(
911 Teuchos::RCP<const std::vector<Real> > vp =
914 std::vector<Real> V(nx_+2,0.0);
915 std::vector<Real> M(nx_+2,0.0);
917 for (
unsigned t = 0; t < nt_+1; t++) {
918 ss = ((t < nt_ && t > 0) ? dt_ : 0.5*dt_);
919 for (
unsigned n = 0; n < nx_+2; n++) {
920 V[n] = (*vp)[t*(nx_+2)+n];
923 for (
unsigned n = 0; n < nx_+2; n++) {
924 (*hvp)[t*(nx_+2)+n] = ss*alpha_*M[n];
Provides the interface to evaluate simulation-based objective functions.
void run_newton(std::vector< Real > &u, const std::vector< Real > &znew, const std::vector< Real > &uold, const std::vector< Real > &zold)
void compute_pde_jacobian(std::vector< Real > &dl, std::vector< Real > &d, std::vector< Real > &du, const std::vector< Real > &u)
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 adjoint of the partial constraint Hessian at , , to vector in direction ...
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)
Real evaluate_target(Real x)
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 adjoint of the partial constraint Hessian at , , to vector in direction ...
void applyInverseJacobian_1(ROL::Vector< Real > &ijv, 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 .
Real dot(const std::vector< Real > &x, const std::vector< Real > &y)
Contains definitions of custom data types in ROL.
void compute_residual(std::vector< Real > &r, const std::vector< Real > &uold, const std::vector< Real > &zold, const std::vector< Real > &unew, const std::vector< Real > &znew)
Real compute_norm(const std::vector< Real > &r)
EqualityConstraint_BurgersControl(int nx=128, int nt=100, Real T=1, Real nu=1.e-2, Real u0=0.0, Real u1=0.0, Real f=0.0)
void solve(ROL::Vector< Real > &c, ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Given , solve for .
virtual void zero()
Set to zero vector.
Defines the linear algebra or vector space interface.
void apply_control_jacobian(std::vector< Real > &jv, const std::vector< Real > &v, bool adjoint=false)
void compute_residual(std::vector< Real > &r, const std::vector< Real > &u, const std::vector< Real > &z)
void apply_pde_hessian(std::vector< Real > &hv, const std::vector< Real > &wold, const std::vector< Real > &vold, const std::vector< Real > &wnew, const std::vector< Real > &vnew)
Defines the equality constraint operator interface for simulation-based optimization.
void hessVec_22(ROL::Vector< Real > &hv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Real value(const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Compute value.
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 applyAdjointJacobian_1(ROL::Vector< Real > &ajv, 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 apply_pde_jacobian_new(std::vector< Real > &jv, const std::vector< Real > &v, const std::vector< Real > &u, bool adjoint=false)
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 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 .
Real dot(const std::vector< Real > &x, const std::vector< Real > &y)
void hessVec_21(ROL::Vector< Real > &hv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
std::vector< Real > u_init_
void hessVec_12(ROL::Vector< Real > &hv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
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.
void scale(std::vector< Real > &u, const Real alpha=0.0)
Objective_BurgersControl(Real alpha=1.e-4, int nx=128, int nt=100, Real T=1.0)
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 .
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 value(ROL::Vector< Real > &c, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Evaluate the constraint operator at .
void applyInverseAdjointJacobian_1(ROL::Vector< Real > &ijv, 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 apply_mass(std::vector< Real > &Mu, const std::vector< Real > &u)
void linear_solve(std::vector< Real > &u, const std::vector< Real > &dl, const std::vector< Real > &d, const std::vector< Real > &du, const std::vector< Real > &r, const bool transpose=false)
void apply_pde_jacobian(std::vector< Real > &jv, const std::vector< Real > &vold, const std::vector< Real > &uold, const std::vector< Real > &vnew, const std::vector< Real > unew, bool adjoint=false)
virtual Real norm() const =0
Returns where .
void update(std::vector< Real > &u, const std::vector< Real > &s, const Real alpha=1.0)
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 adjoint of the partial constraint Hessian at , , to vector in direction ...
void apply_pde_jacobian_old(std::vector< Real > &jv, const std::vector< Real > &v, const std::vector< Real > &u, bool adjoint=false)
static const double ROL_EPSILON
Platform-dependent machine epsilon.
void applyAdjointHessian_11(ROL::Vector< Real > &hwv, const ROL::Vector< Real > &w, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Apply the adjoint of the partial constraint Hessian at , , to vector in direction ...