59 #include "ROL_ParameterList.hpp"
62 #include "Teuchos_GlobalMPISession.hpp"
63 #include "Teuchos_LAPACK.hpp"
87 return std::sqrt(
dot(r,r));
90 Real
dot(
const std::vector<Real> &x,
const std::vector<Real> &y) {
92 Real c = ((x.size()==
nx_) ? 4.0 : 2.0);
93 for (
unsigned i = 0; i < x.size(); i++) {
95 ip +=
dx_/6.0*(c*x[i] + x[i+1])*y[i];
97 else if ( i == x.size()-1 ) {
98 ip +=
dx_/6.0*(x[i-1] + c*x[i])*y[i];
101 ip +=
dx_/6.0*(x[i-1] + 4.0*x[i] + x[i+1])*y[i];
109 void update(std::vector<Real> &u,
const std::vector<Real> &s,
const Real alpha=1.0) {
110 for (
unsigned i = 0; i < u.size(); i++) {
115 void scale(std::vector<Real> &u,
const Real alpha=0.0) {
116 for (
unsigned i = 0; i < u.size(); i++) {
121 void compute_residual(std::vector<Real> &r,
const std::vector<Real> &uold,
const std::vector<Real> &zold,
122 const std::vector<Real> &unew,
const std::vector<Real> &znew) {
125 for (
unsigned n = 0; n <
nx_; n++) {
139 r[n] -= 0.5*
dt_*unew[n-1]*(unew[n-1]+unew[n])/6.0;
140 r[n] -= 0.5*
dt_*uold[n-1]*(uold[n-1]+uold[n])/6.0;
143 r[n] += 0.5*
dt_*unew[n+1]*(unew[n]+unew[n+1])/6.0;
144 r[n] += 0.5*
dt_*uold[n+1]*(uold[n]+uold[n+1])/6.0;
147 r[n] -= 0.5*
dt_*
dx_/6.0*(znew[n]+4.0*znew[n+1]+znew[n+2]);
148 r[n] -= 0.5*
dt_*
dx_/6.0*(zold[n]+4.0*zold[n+1]+zold[n+2]);
158 const std::vector<Real> &u) {
167 for (
unsigned n = 0; n <
nx_; n++) {
169 dl[n] += 0.5*
dt_*(-2.0*u[n]-u[n+1])/6.0;
170 d[n] += 0.5*
dt_*u[n+1]/6.0;
173 d[n] -= 0.5*
dt_*u[n-1]/6.0;
174 du[n-1] += 0.5*
dt_*(u[n-1]+2.0*u[n])/6.0;
179 d[nx_-1] += 0.5*
dt_*
u1_/6.0;
183 bool adjoint =
false) {
187 for (
unsigned n = 0; n <
nx_; n++) {
192 jv[n] -= 0.5*
dt_*(u[n-1]/6.0*v[n]-(u[n-1]+2.0*u[n])/6.0*v[n-1]);
195 jv[n] -= 0.5*
dt_*(u[n-1]/6.0*v[n]+(u[n]+2.0*u[n-1])/6.0*v[n-1]);
201 jv[n] += 0.5*
dt_*(u[n+1]/6.0*v[n]-(u[n+1]+2.0*u[n])/6.0*v[n+1]);
204 jv[n] += 0.5*
dt_*(u[n+1]/6.0*v[n]+(u[n]+2.0*u[n+1])/6.0*v[n+1]);
208 jv[0] -= 0.5*
dt_*
u0_/6.0*v[0];
209 jv[nx_-1] += 0.5*
dt_*
u1_/6.0*v[nx_-1];
213 bool adjoint =
false) {
217 for (
unsigned n = 0; n <
nx_; n++) {
222 jv[n] -= 0.5*
dt_*(u[n-1]/6.0*v[n]-(u[n-1]+2.0*u[n])/6.0*v[n-1]);
225 jv[n] -= 0.5*
dt_*(u[n-1]/6.0*v[n]+(u[n]+2.0*u[n-1])/6.0*v[n-1]);
231 jv[n] += 0.5*
dt_*(u[n+1]/6.0*v[n]-(u[n+1]+2.0*u[n])/6.0*v[n+1]);
234 jv[n] += 0.5*
dt_*(u[n+1]/6.0*v[n]+(u[n]+2.0*u[n+1])/6.0*v[n+1]);
238 jv[0] -= 0.5*
dt_*
u0_/6.0*v[0];
239 jv[nx_-1] += 0.5*
dt_*
u1_/6.0*v[nx_-1];
242 void apply_pde_jacobian(std::vector<Real> &jv,
const std::vector<Real> &vold,
const std::vector<Real> &uold,
243 const std::vector<Real> &vnew,
const std::vector<Real> unew,
bool adjoint =
false) {
247 for (
unsigned n = 0; n <
nx_; n++) {
254 jv[n] -= 0.5*
dt_*(unew[n-1]/6.0*vnew[n]-(unew[n-1]+2.0*unew[n])/6.0*vnew[n-1]);
255 jv[n] -= 0.5*
dt_*(uold[n-1]/6.0*vold[n]-(uold[n-1]+2.0*uold[n])/6.0*vold[n-1]);
258 jv[n] -= 0.5*
dt_*(unew[n-1]/6.0*vnew[n]+(unew[n]+2.0*unew[n-1])/6.0*vnew[n-1]);
259 jv[n] -= 0.5*
dt_*(uold[n-1]/6.0*vold[n]+(uold[n]+2.0*uold[n-1])/6.0*vold[n-1]);
266 jv[n] += 0.5*
dt_*(unew[n+1]/6.0*vnew[n]-(unew[n+1]+2.0*unew[n])/6.0*vnew[n+1]);
267 jv[n] += 0.5*
dt_*(uold[n+1]/6.0*vold[n]-(uold[n+1]+2.0*uold[n])/6.0*vold[n+1]);
270 jv[n] += 0.5*
dt_*(unew[n+1]/6.0*vnew[n]+(unew[n]+2.0*unew[n+1])/6.0*vnew[n+1]);
271 jv[n] += 0.5*
dt_*(uold[n+1]/6.0*vold[n]+(uold[n]+2.0*uold[n+1])/6.0*vold[n+1]);
275 jv[0] -= 0.5*
dt_*
u0_/6.0*vnew[0];
276 jv[0] -= 0.5*
dt_*
u0_/6.0*vold[0];
277 jv[nx_-1] += 0.5*
dt_*
u1_/6.0*vnew[nx_-1];
278 jv[nx_-1] += 0.5*
dt_*
u1_/6.0*vold[nx_-1];
281 void apply_pde_hessian(std::vector<Real> &hv,
const std::vector<Real> &wold,
const std::vector<Real> &vold,
282 const std::vector<Real> &wnew,
const std::vector<Real> &vnew) {
285 for (
unsigned n = 0; n <
nx_; n++) {
287 hv[n] += 0.5*
dt_*((wnew[n-1]*(vnew[n-1]+2.0*vnew[n]) - wnew[n]*vnew[n-1])/6.0);
288 hv[n] += 0.5*
dt_*((wold[n-1]*(vold[n-1]+2.0*vold[n]) - wold[n]*vold[n-1])/6.0);
291 hv[n] += 0.5*
dt_*((wnew[n]*vnew[n+1] - wnew[n+1]*(2.0*vnew[n]+vnew[n+1]))/6.0);
292 hv[n] += 0.5*
dt_*((wold[n]*vold[n+1] - wold[n+1]*(2.0*vold[n]+vold[n+1]))/6.0);
299 unsigned dim = ((adjoint ==
true) ?
nx_+2 :
nx_);
301 for (
unsigned n = 0; n <
dim; n++) {
304 jv[n] = -0.5*
dt_*
dx_/6.0*v[n];
307 jv[n] = -0.5*
dt_*
dx_/6.0*(4.0*v[n-1]+v[n]);
309 else if ( n ==
nx_ ) {
310 jv[n] = -0.5*
dt_*
dx_/6.0*(4.0*v[n-1]+v[n-2]);
312 else if ( n ==
nx_+1 ) {
313 jv[n] = -0.5*
dt_*
dx_/6.0*v[n-2];
316 jv[n] = -0.5*
dt_*
dx_/6.0*(v[n-2]+4.0*v[n-1]+v[n]);
320 jv[n] -= 0.5*
dt_*
dx_/6.0*(v[n]+4.0*v[n+1]+v[n+2]);
325 void run_newton(std::vector<Real> &u,
const std::vector<Real> &znew,
326 const std::vector<Real> &uold,
const std::vector<Real> &zold) {
330 std::vector<Real> r(
nx_,0.0);
334 Real rtol = 1.e2*ROL::ROL_EPSILON<Real>();
335 unsigned maxit = 500;
337 std::vector<Real> d(
nx_,0.0);
338 std::vector<Real> dl(
nx_-1,0.0);
339 std::vector<Real> du(
nx_-1,0.0);
341 Real alpha = 1.0, tmp = 0.0;
342 std::vector<Real> s(
nx_,0.0);
343 std::vector<Real> utmp(
nx_,0.0);
344 for (
unsigned i = 0; i < maxit; i++) {
353 utmp.assign(u.begin(),u.end());
357 while ( rnorm > (1.0-1.e-4*alpha)*tmp && alpha > std::sqrt(ROL::ROL_EPSILON<Real>()) ) {
359 utmp.assign(u.begin(),u.end());
365 u.assign(utmp.begin(),utmp.end());
366 if ( rnorm < rtol ) {
373 const std::vector<Real> &dl,
const std::vector<Real> &d,
const std::vector<Real> &du,
374 const std::vector<Real> &r,
const bool transpose =
false) {
375 bool useLAPACK =
true;
377 u.assign(r.begin(),r.end());
379 std::vector<Real> Dl(dl);
380 std::vector<Real> Du(du);
381 std::vector<Real> D(d);
383 Teuchos::LAPACK<int,Real> lp;
384 std::vector<Real> Du2(
nx_-2,0.0);
385 std::vector<int> ipiv(
nx_,0);
389 lp.GTTRF(
nx_,&Dl[0],&D[0],&Du[0],&Du2[0],&ipiv[0],&info);
390 char trans = ((transpose ==
true) ?
'T' :
'N');
391 lp.GTTRS(trans,
nx_,nhrs,&Dl[0],&D[0],&Du[0],&Du2[0],&ipiv[0],&u[0],ldb,&info);
396 unsigned maxit = 100;
397 Real rtol = std::min(1.e-12,1.e-4*std::sqrt(
dot(r,r)));
400 Real rnorm = 10.0*rtol;
401 for (
unsigned i = 0; i < maxit; i++) {
402 for (
unsigned n = 0; n <
nx_; n++) {
405 u[n] -= ((transpose ==
false) ? du[n] : dl[n])*u[n+1]/d[n];
408 u[n] -= ((transpose ==
false) ? dl[n-1] : du[n-1])*u[n-1]/d[n];
413 for (
unsigned n = 0; n <
nx_; n++) {
414 resid = r[n] - d[n]*u[n];
416 resid -= ((transpose ==
false) ? du[n] : dl[n])*u[n+1];
419 resid -= ((transpose ==
false) ? dl[n-1] : du[n-1])*u[n-1];
421 rnorm += resid*resid;
423 rnorm = std::sqrt(rnorm);
424 if ( rnorm < rtol ) {
435 Real nu = 1.e-2, Real u0 = 0.0, Real u1 = 0.0, Real f = 0.0)
437 dx_ = 1.0/((Real)nx+1.0);
442 for (
unsigned n = 0; n <
nx_; n++) {
444 u_init_[n] = ((x <= 0.5) ? 1.0 : 0.0);
449 ROL::Ptr<std::vector<Real> > cp =
451 ROL::Ptr<const std::vector<Real> > up =
453 ROL::Ptr<const std::vector<Real> > zp =
456 std::vector<Real> C(
nx_,0.0);
457 std::vector<Real> uold(
u_init_);
458 std::vector<Real> unew(
u_init_);
459 std::vector<Real> zold(
nx_+2,0.0);
460 std::vector<Real> znew(
nx_+2,0.0);
462 for (
unsigned n = 0; n <
nx_+2; n++) {
466 for (
unsigned t = 0; t <
nt_; t++) {
468 for (
unsigned n = 0; n <
nx_; n++) {
469 unew[n] = (*up)[t*nx_+n];
472 for (
unsigned n = 0; n < nx_+2; n++) {
473 znew[n] = (*zp)[(t+1)*(nx_+2)+n];
478 for (
unsigned n = 0; n <
nx_; n++) {
479 (*cp)[t*nx_+n] = C[n];
482 uold.assign(unew.begin(),unew.end());
483 zold.assign(znew.begin(),znew.end());
488 ROL::Ptr<std::vector<Real> > up =
490 up->assign(up->size(),z.
norm()/up->size());
491 ROL::Ptr<const std::vector<Real> > zp =
494 std::vector<Real> uold(
u_init_);
495 std::vector<Real> unew(
u_init_);
496 std::vector<Real> zold(
nx_+2,0.0);
497 std::vector<Real> znew(
nx_+2,0.0);
499 for (
unsigned n = 0; n <
nx_+2; n++) {
503 for (
unsigned t = 0; t <
nt_; t++) {
505 for (
unsigned n = 0; n < nx_+2; n++) {
506 znew[n] = (*zp)[(t+1)*(nx_+2)+n];
511 for (
unsigned n = 0; n <
nx_; n++) {
512 (*up)[t*nx_+n] = unew[n];
515 uold.assign(unew.begin(),unew.end());
516 zold.assign(znew.begin(),znew.end());
523 ROL::Ptr<std::vector<Real> > jvp =
525 ROL::Ptr<const std::vector<Real> > vp =
527 ROL::Ptr<const std::vector<Real> > up =
529 std::vector<Real> J(
nx_,0.0);
530 std::vector<Real> vold(
nx_,0.0);
531 std::vector<Real> vnew(
nx_,0.0);
532 std::vector<Real> uold(
nx_,0.0);
533 std::vector<Real> unew(
nx_,0.0);
534 for (
unsigned t = 0; t <
nt_; t++) {
535 for (
unsigned n = 0; n <
nx_; n++) {
536 unew[n] = (*up)[t*nx_+n];
537 vnew[n] = (*vp)[t*nx_+n];
540 for (
unsigned n = 0; n <
nx_; n++) {
541 (*jvp)[t*nx_+n] = J[n];
543 vold.assign(vnew.begin(),vnew.end());
544 uold.assign(unew.begin(),unew.end());
550 ROL::Ptr<std::vector<Real> > jvp =
552 ROL::Ptr<const std::vector<Real> > vp =
554 ROL::Ptr<const std::vector<Real> > zp =
556 std::vector<Real> vnew(
nx_+2,0.0);
557 std::vector<Real> vold(
nx_+2,0.0);
558 std::vector<Real> jnew(
nx_,0.0);
559 std::vector<Real> jold(
nx_,0.0);
560 for (
unsigned n = 0; n <
nx_+2; n++) {
564 for (
unsigned t = 0; t <
nt_; t++) {
565 for (
unsigned n = 0; n < nx_+2; n++) {
566 vnew[n] = (*vp)[(t+1)*(nx_+2)+n];
569 for (
unsigned n = 0; n <
nx_; n++) {
571 (*jvp)[t*nx_+n] = jnew[n] + jold[n];
573 jold.assign(jnew.begin(),jnew.end());
579 ROL::Ptr<std::vector<Real> > ijvp =
581 ROL::Ptr<const std::vector<Real> > vp =
583 ROL::Ptr<const std::vector<Real> > up =
585 std::vector<Real> J(
nx_,0.0);
586 std::vector<Real> r(
nx_,0.0);
587 std::vector<Real> s(
nx_,0.0);
588 std::vector<Real> uold(
nx_,0.0);
589 std::vector<Real> unew(
nx_,0.0);
590 std::vector<Real> d(
nx_,0.0);
591 std::vector<Real> dl(
nx_-1,0.0);
592 std::vector<Real> du(
nx_-1,0.0);
593 for (
unsigned t = 0; t <
nt_; t++) {
595 for (
unsigned n = 0; n <
nx_; n++) {
596 r[n] = (*vp)[t*nx_+n];
597 unew[n] = (*up)[t*nx_+n];
606 for (
unsigned n = 0; n <
nx_; n++) {
607 (*ijvp)[t*nx_+n] = s[n];
609 uold.assign(unew.begin(),unew.end());
615 ROL::Ptr<std::vector<Real> > jvp =
617 ROL::Ptr<const std::vector<Real> > vp =
619 ROL::Ptr<const std::vector<Real> > up =
621 std::vector<Real> J(
nx_,0.0);
622 std::vector<Real> vold(
nx_,0.0);
623 std::vector<Real> vnew(
nx_,0.0);
624 std::vector<Real> unew(
nx_,0.0);
625 for (
unsigned t =
nt_; t > 0; t--) {
626 for (
unsigned n = 0; n <
nx_; n++) {
627 unew[n] = (*up)[(t-1)*nx_+n];
628 vnew[n] = (*vp)[(t-1)*nx_+n];
631 for (
unsigned n = 0; n <
nx_; n++) {
632 (*jvp)[(t-1)*nx_+n] = J[n];
634 vold.assign(vnew.begin(),vnew.end());
640 ROL::Ptr<std::vector<Real> > jvp =
642 ROL::Ptr<const std::vector<Real> > vp =
644 ROL::Ptr<const std::vector<Real> > zp =
646 std::vector<Real> vnew(
nx_,0.0);
647 std::vector<Real> vold(
nx_,0.0);
648 std::vector<Real> jnew(
nx_+2,0.0);
649 std::vector<Real> jold(
nx_+2,0.0);
650 for (
unsigned t =
nt_+1; t > 0; t--) {
651 for (
unsigned n = 0; n <
nx_; n++) {
653 vnew[n] = (*vp)[(t-2)*nx_+n];
660 for (
unsigned n = 0; n < nx_+2; n++) {
662 (*jvp)[(t-1)*(nx_+2)+n] = jnew[n] + jold[n];
664 jold.assign(jnew.begin(),jnew.end());
670 ROL::Ptr<std::vector<Real> > ijvp =
672 ROL::Ptr<const std::vector<Real> > vp =
674 ROL::Ptr<const std::vector<Real> > up =
676 std::vector<Real> J(
nx_,0.0);
677 std::vector<Real> r(
nx_,0.0);
678 std::vector<Real> s(
nx_,0.0);
679 std::vector<Real> unew(
nx_,0.0);
680 std::vector<Real> d(
nx_,0.0);
681 std::vector<Real> dl(
nx_-1,0.0);
682 std::vector<Real> du(
nx_-1,0.0);
683 for (
unsigned t =
nt_; t > 0; t--) {
685 for (
unsigned n = 0; n <
nx_; n++) {
686 r[n] = (*vp)[(t-1)*nx_+n];
687 unew[n] = (*up)[(t-1)*nx_+n];
696 for (
unsigned n = 0; n <
nx_; n++) {
697 (*ijvp)[(t-1)*nx_+n] = s[n];
704 ROL::Ptr<std::vector<Real> > hwvp =
706 ROL::Ptr<const std::vector<Real> > wp =
708 ROL::Ptr<const std::vector<Real> > vp =
710 std::vector<Real> snew(
nx_,0.0);
711 std::vector<Real> wnew(
nx_,0.0);
712 std::vector<Real> wold(
nx_,0.0);
713 std::vector<Real> vnew(
nx_,0.0);
714 for (
unsigned t =
nt_; t > 0; t--) {
715 for (
unsigned n = 0; n <
nx_; n++) {
716 vnew[n] = (*vp)[(t-1)*nx_+n];
717 wnew[n] = (*wp)[(t-1)*nx_+n];
720 for (
unsigned n = 0; n <
nx_; n++) {
721 (*hwvp)[(t-1)*nx_+n] = snew[n];
723 wold.assign(wnew.begin(),wnew.end());
759 case 1: val = ((x<=0.5) ? 1.0 : 0.0);
break;
760 case 2: val = 1.0;
break;
761 case 3: val = std::abs(std::sin(8.0*M_PI*x));
break;
762 case 4: val = std::exp(-0.5*(x-0.5)*(x-0.5));
break;
767 Real
dot(
const std::vector<Real> &x,
const std::vector<Real> &y) {
769 Real c = ((x.size()==
nx_) ? 4.0 : 2.0);
770 for (
unsigned i=0; i<x.size(); i++) {
772 ip +=
dx_/6.0*(c*x[i] + x[i+1])*y[i];
774 else if ( i == x.size()-1 ) {
775 ip +=
dx_/6.0*(x[i-1] + c*x[i])*y[i];
778 ip +=
dx_/6.0*(x[i-1] + 4.0*x[i] + x[i+1])*y[i];
784 void apply_mass(std::vector<Real> &Mu,
const std::vector<Real> &u ) {
785 Mu.resize(u.size(),0.0);
786 Real c = ((u.size()==
nx_) ? 4.0 : 2.0);
787 for (
unsigned i=0; i<u.size(); i++) {
789 Mu[i] =
dx_/6.0*(c*u[i] + u[i+1]);
791 else if ( i == u.size()-1 ) {
792 Mu[i] =
dx_/6.0*(u[i-1] + c*u[i]);
795 Mu[i] =
dx_/6.0*(u[i-1] + 4.0*u[i] + u[i+1]);
807 dx_ = 1.0/((Real)nx+1.0);
812 ROL::Ptr<const std::vector<Real> > up =
814 ROL::Ptr<const std::vector<Real> > zp =
817 std::vector<Real> U(
nx_,0.0);
818 std::vector<Real> G(
nx_,0.0);
819 std::vector<Real> Z(
nx_+2,0.0);
820 for (
unsigned n = 0; n <
nx_+2; n++) {
825 for (
unsigned t = 0; t <
nt_; t++) {
826 ss = ((t < nt_-1) ?
dt_ : 0.5*
dt_);
827 for (
unsigned n = 0; n <
nx_; n++) {
831 val += 0.5*ss*
dot(U,U);
832 val -= 0.5*ss*
dot(G,G);
833 for (
unsigned n = 0; n < nx_+2; n++) {
834 Z[n] = (*zp)[(t+1)*(nx_+2)+n];
842 ROL::Ptr<std::vector<Real> > gp =
844 ROL::Ptr<const std::vector<Real> > up =
847 std::vector<Real> U(
nx_,0.0);
848 std::vector<Real> M(
nx_,0.0);
850 for (
unsigned t = 0; t <
nt_; t++) {
851 ss = ((t < nt_-1) ?
dt_ : 0.5*
dt_);
852 for (
unsigned n = 0; n <
nx_; n++) {
856 for (
unsigned n = 0; n <
nx_; n++) {
857 (*gp)[t*nx_+n] = ss*M[n];
863 ROL::Ptr<std::vector<Real> > gp =
865 ROL::Ptr<const std::vector<Real> > zp =
868 std::vector<Real> Z(
nx_+2,0.0);
869 std::vector<Real> M(
nx_+2,0.0);
871 for (
unsigned t = 0; t <
nt_+1; t++) {
872 ss = ((t < nt_ && t > 0) ?
dt_ : 0.5*
dt_);
873 for (
unsigned n = 0; n <
nx_+2; n++) {
874 Z[n] = (*zp)[t*(nx_+2)+n];
877 for (
unsigned n = 0; n < nx_+2; n++) {
878 (*gp)[t*(nx_+2)+n] = ss*
alpha_*M[n];
885 ROL::Ptr<std::vector<Real> > hvp =
887 ROL::Ptr<const std::vector<Real> > vp =
890 std::vector<Real>
V(
nx_,0.0);
891 std::vector<Real> M(
nx_,0.0);
893 for (
unsigned t = 0; t <
nt_; t++) {
894 ss = ((t < nt_-1) ?
dt_ : 0.5*
dt_);
895 for (
unsigned n = 0; n <
nx_; n++) {
896 V[n] = (*vp)[t*nx_+n];
899 for (
unsigned n = 0; n <
nx_; n++) {
900 (*hvp)[t*nx_+n] = ss*M[n];
917 ROL::Ptr<std::vector<Real> > hvp = ROL::constPtrCast<std::vector<Real> >(
919 ROL::Ptr<const std::vector<Real> > vp =
922 std::vector<Real>
V(
nx_+2,0.0);
923 std::vector<Real> M(
nx_+2,0.0);
925 for (
unsigned t = 0; t <
nt_+1; t++) {
926 ss = ((t < nt_ && t > 0) ?
dt_ : 0.5*
dt_);
927 for (
unsigned n = 0; n <
nx_+2; n++) {
928 V[n] = (*vp)[t*(nx_+2)+n];
931 for (
unsigned n = 0; n < nx_+2; n++) {
932 (*hvp)[t*(nx_+2)+n] = ss*
alpha_*M[n];
void apply_pde_jacobian_old(std::vector< Real > &jv, const std::vector< Real > &v, const std::vector< Real > &u, bool adjoint=false)
Provides the interface to evaluate simulation-based objective functions.
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...
Real evaluate_target(Real x)
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...
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 simulation-space derivative of the adjoint of the constraint simulation-space Jacobian at ...
void run_newton(std::vector< Real > &u, const std::vector< Real > &znew, const std::vector< Real > &uold, const std::vector< Real > &zold)
void solve(ROL::Vector< Real > &c, ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Given , solve for .
Contains definitions of custom data types in ROL.
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 ...
virtual void zero()
Set to zero vector.
Defines the linear algebra or vector space interface.
Defines a no-output stream class ROL::NullStream and a function makeStreamPtr which either wraps a re...
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 apply_pde_jacobian_new(std::vector< Real > &jv, const std::vector< Real > &v, const std::vector< Real > &u, bool adjoint=false)
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 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)
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 compute_pde_jacobian(std::vector< Real > &dl, std::vector< Real > &d, std::vector< Real > &du, const std::vector< Real > &u)
Real dot(const std::vector< Real > &x, const std::vector< Real > &y)
Constraint_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 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 hessVec_21(ROL::Vector< Real > &hv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
void apply_control_jacobian(std::vector< Real > &jv, const std::vector< Real > &v, bool adjoint=false)
void hessVec_12(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 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.
Real dot(const std::vector< Real > &x, const std::vector< Real > &y)
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 .
void compute_residual(std::vector< Real > &r, const std::vector< Real > &u, const std::vector< Real > &z)
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 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 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 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_mass(std::vector< Real > &Mu, const std::vector< Real > &u)
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 ...
virtual Real norm() const =0
Returns where .
Real compute_norm(const std::vector< Real > &r)
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)
void update(std::vector< Real > &u, const std::vector< Real > &s, const Real alpha=1.0)
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...
Defines the constraint operator interface for simulation-based optimization.
void value(ROL::Vector< Real > &c, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Evaluate the constraint operator at .
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 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)