55 #include "Teuchos_LAPACK.hpp" 83 void update(std::vector<Real> &u,
const std::vector<Real> &s,
const Real alpha=1.0)
const {
84 for (
unsigned i=0; i<u.size(); i++) {
89 void axpy(std::vector<Real> &out,
const Real a,
const std::vector<Real> &x,
const std::vector<Real> &y)
const {
90 for (
unsigned i=0; i < x.size(); i++) {
91 out[i] = a*x[i] + y[i];
95 void scale(std::vector<Real> &u,
const Real alpha=0.0)
const {
96 for (
unsigned i=0; i<u.size(); i++) {
101 void linear_solve(std::vector<Real> &u, std::vector<Real> &dl, std::vector<Real> &d, std::vector<Real> &du,
102 const std::vector<Real> &r,
const bool transpose =
false)
const {
103 if ( r.size() == 1 ) {
104 u.resize(1,r[0]/d[0]);
106 else if ( r.size() == 2 ) {
108 Real det = d[0]*d[1] - dl[0]*du[0];
109 u[0] = (d[1]*r[0] - du[0]*r[1])/det;
110 u[1] = (d[0]*r[1] - dl[0]*r[0])/det;
113 u.assign(r.begin(),r.end());
115 Teuchos::LAPACK<int,Real> lp;
116 std::vector<Real> du2(r.size()-2,0.0);
117 std::vector<int> ipiv(r.size(),0);
122 lp.GTTRF(dim,&dl[0],&d[0],&du[0],&du2[0],&ipiv[0],&info);
127 lp.GTTRS(trans,dim,nhrs,&dl[0],&d[0],&du[0],&du2[0],&ipiv[0],&u[0],ldb,&info);
132 BurgersFEM(
int nx = 128, Real nl = 1.0, Real cH1 = 1.0, Real cL2 = 1.0)
133 : nx_(nx), dx_(1.0/((Real)nx+1.0)), nl_(nl), cH1_(cH1), cL2_(cL2) {
141 nu_ = std::pow(10.0,nu-2.0);
159 Real
compute_L2_dot(
const std::vector<Real> &x,
const std::vector<Real> &y)
const {
161 Real c = (((int)x.size()==
nx_) ? 4.0 : 2.0);
162 for (
unsigned i=0; i<x.size(); i++) {
164 ip += dx_/6.0*(c*x[i] + x[i+1])*y[i];
166 else if ( i == x.size()-1 ) {
167 ip += dx_/6.0*(x[i-1] + c*x[i])*y[i];
170 ip += dx_/6.0*(x[i-1] + 4.0*x[i] + x[i+1])*y[i];
182 void apply_mass(std::vector<Real> &Mu,
const std::vector<Real> &u )
const {
183 Mu.resize(u.size(),0.0);
184 Real c = (((int)u.size()==
nx_) ? 4.0 : 2.0);
185 for (
unsigned i=0; i<u.size(); i++) {
187 Mu[i] = dx_/6.0*(c*u[i] + u[i+1]);
189 else if ( i == u.size()-1 ) {
190 Mu[i] = dx_/6.0*(u[i-1] + c*u[i]);
193 Mu[i] = dx_/6.0*(u[i-1] + 4.0*u[i] + u[i+1]);
200 unsigned nx = u.size();
202 std::vector<Real> dl(nx-1,dx_/6.0);
203 std::vector<Real> d(nx,2.0*dx_/3.0);
204 std::vector<Real> du(nx-1,dx_/6.0);
205 if ( (
int)nx != nx_ ) {
214 std::vector<Real> u(nx_,0.0), Mu(nx_,0.0), iMMu(nx_,0.0), diff(nx_,0.0);
215 for (
int i = 0; i <
nx_; i++) {
216 u[i] = 2.0*(Real)rand()/(Real)RAND_MAX - 1.0;
220 axpy(diff,-1.0,iMMu,u);
223 outStream <<
"\nTest Inverse State Mass Matrix\n";
224 outStream <<
" ||u - inv(M)Mu|| = " << error <<
"\n";
225 outStream <<
" ||u|| = " << normu <<
"\n";
226 outStream <<
" Relative Error = " << error/normu <<
"\n";
229 u.resize(nx_+2,0.0); Mu.resize(nx_+2,0.0); iMMu.resize(nx_+2,0.0); diff.resize(nx_+2,0.0);
230 for (
int i = 0; i < nx_+2; i++) {
231 u[i] = 2.0*(Real)rand()/(Real)RAND_MAX - 1.0;
235 axpy(diff,-1.0,iMMu,u);
238 outStream <<
"\nTest Inverse Control Mass Matrix\n";
239 outStream <<
" ||z - inv(M)Mz|| = " << error <<
"\n";
240 outStream <<
" ||z|| = " << normu <<
"\n";
241 outStream <<
" Relative Error = " << error/normu <<
"\n";
249 Real
compute_H1_dot(
const std::vector<Real> &x,
const std::vector<Real> &y)
const {
251 for (
int i=0; i<
nx_; i++) {
253 ip += cL2_*dx_/6.0*(4.0*x[i] + x[i+1])*y[i];
254 ip += cH1_*(2.0*x[i] - x[i+1])/dx_*y[i];
256 else if ( i == nx_-1 ) {
257 ip += cL2_*dx_/6.0*(x[i-1] + 4.0*x[i])*y[i];
258 ip += cH1_*(2.0*x[i] - x[i-1])/dx_*y[i];
261 ip += cL2_*dx_/6.0*(x[i-1] + 4.0*x[i] + x[i+1])*y[i];
262 ip += cH1_*(2.0*x[i] - x[i-1] - x[i+1])/dx_*y[i];
274 void apply_H1(std::vector<Real> &Mu,
const std::vector<Real> &u )
const {
276 for (
int i=0; i<
nx_; i++) {
278 Mu[i] = cL2_*dx_/6.0*(4.0*u[i] + u[i+1])
279 + cH1_*(2.0*u[i] - u[i+1])/
dx_;
281 else if ( i == nx_-1 ) {
282 Mu[i] = cL2_*dx_/6.0*(u[i-1] + 4.0*u[i])
283 + cH1_*(2.0*u[i] - u[i-1])/
dx_;
286 Mu[i] = cL2_*dx_/6.0*(u[i-1] + 4.0*u[i] + u[i+1])
287 + cH1_*(2.0*u[i] - u[i-1] - u[i+1])/
dx_;
295 std::vector<Real> dl(nx_-1,cL2_*dx_/6.0 - cH1_/dx_);
296 std::vector<Real> d(nx_,2.0*(cL2_*dx_/3.0 + cH1_/dx_));
297 std::vector<Real> du(nx_-1,cL2_*dx_/6.0 - cH1_/dx_);
302 std::vector<Real> u(nx_,0.0), Mu(nx_,0.0), iMMu(nx_,0.0), diff(nx_,0.0);
303 for (
int i = 0; i <
nx_; i++) {
304 u[i] = 2.0*(Real)rand()/(Real)RAND_MAX - 1.0;
308 axpy(diff,-1.0,iMMu,u);
311 outStream <<
"\nTest Inverse State H1 Matrix\n";
312 outStream <<
" ||u - inv(M)Mu|| = " << error <<
"\n";
313 outStream <<
" ||u|| = " << normu <<
"\n";
314 outStream <<
" Relative Error = " << error/normu <<
"\n";
323 const std::vector<Real> &z)
const {
326 for (
int i=0; i<
nx_; i++) {
329 r[i] = nu_/dx_*(2.0*u[i]-u[i+1]);
332 r[i] = nu_/dx_*(2.0*u[i]-u[i-1]);
335 r[i] = nu_/dx_*(2.0*u[i]-u[i-1]-u[i+1]);
339 r[i] += nl_*u[i+1]*(u[i]+u[i+1])/6.0;
342 r[i] -= nl_*u[i-1]*(u[i-1]+u[i])/6.0;
345 r[i] -= dx_/6.0*(z[i]+4.0*z[i+1]+z[i+2]);
350 r[0] -= nl_*(u0_*u[ 0]/6.0 + u0_*u0_/6.0) + nu_*u0_/
dx_;
351 r[nx_-1] += nl_*(u1_*u[nx_-1]/6.0 + u1_*u1_/6.0) - nu_*u1_/
dx_;
359 std::vector<Real> &d,
360 std::vector<Real> &du,
361 const std::vector<Real> &u)
const {
364 d.resize(nx_,nu_*2.0/dx_);
366 dl.resize(nx_-1,-nu_/dx_);
368 du.resize(nx_-1,-nu_/dx_);
370 for (
int i=0; i<
nx_; i++) {
372 dl[i] += nl_*(-2.0*u[i]-u[i+1])/6.0;
373 d[i] += nl_*u[i+1]/6.0;
376 d[i] -= nl_*u[i-1]/6.0;
377 du[i-1] += nl_*(u[i-1]+2.0*u[i])/6.0;
381 d[ 0] -= nl_*u0_/6.0;
382 d[nx_-1] += nl_*u1_/6.0;
387 const std::vector<Real> &v,
388 const std::vector<Real> &u,
389 const std::vector<Real> &z)
const {
391 for (
int i = 0; i <
nx_; i++) {
392 jv[i] = nu_/dx_*2.0*v[i];
394 jv[i] += -nu_/dx_*v[i-1]-nl_*(u[i-1]/6.0*v[i]+(u[i]+2.0*u[i-1])/6.0*v[i-1]);
397 jv[i] += -nu_/dx_*v[i+1]+nl_*(u[i+1]/6.0*v[i]+(u[i]+2.0*u[i+1])/6.0*v[i+1]);
400 jv[ 0] -= nl_*u0_/6.0*v[0];
401 jv[nx_-1] += nl_*u1_/6.0*v[nx_-1];
406 const std::vector<Real> &v,
407 const std::vector<Real> &u,
408 const std::vector<Real> &z)
const {
410 std::vector<Real> d(nx_,0.0);
411 std::vector<Real> dl(nx_-1,0.0);
412 std::vector<Real> du(nx_-1,0.0);
420 const std::vector<Real> &v,
421 const std::vector<Real> &u,
422 const std::vector<Real> &z)
const {
424 for (
int i = 0; i <
nx_; i++) {
425 ajv[i] = nu_/dx_*2.0*v[i];
427 ajv[i] += -nu_/dx_*v[i-1]-nl_*(u[i-1]/6.0*v[i]
428 -(u[i-1]+2.0*u[i])/6.0*v[i-1]);
431 ajv[i] += -nu_/dx_*v[i+1]+nl_*(u[i+1]/6.0*v[i]
432 -(u[i+1]+2.0*u[i])/6.0*v[i+1]);
435 ajv[ 0] -= nl_*u0_/6.0*v[0];
436 ajv[nx_-1] += nl_*u1_/6.0*v[nx_-1];
441 const std::vector<Real> &v,
442 const std::vector<Real> &u,
443 const std::vector<Real> &z)
const {
445 std::vector<Real> d(nx_,0.0);
446 std::vector<Real> du(nx_-1,0.0);
447 std::vector<Real> dl(nx_-1,0.0);
458 const std::vector<Real> &v,
459 const std::vector<Real> &u,
460 const std::vector<Real> &z)
const {
461 for (
int i=0; i<
nx_; i++) {
463 jv[i] = -dx_/6.0*(v[i]+4.0*v[i+1]+v[i+2]);
469 const std::vector<Real> &v,
470 const std::vector<Real> &u,
471 const std::vector<Real> &z)
const {
472 for (
int i=0; i<nx_+2; i++) {
474 jv[i] = -dx_/6.0*v[i];
477 jv[i] = -dx_/6.0*(4.0*v[i-1]+v[i]);
479 else if ( i == nx_ ) {
480 jv[i] = -dx_/6.0*(4.0*v[i-1]+v[i-2]);
482 else if ( i == nx_+1 ) {
483 jv[i] = -dx_/6.0*v[i-2];
486 jv[i] = -dx_/6.0*(v[i-2]+4.0*v[i-1]+v[i]);
495 const std::vector<Real> &w,
496 const std::vector<Real> &v,
497 const std::vector<Real> &u,
498 const std::vector<Real> &z)
const {
499 for (
int i=0; i<
nx_; i++) {
503 ahwv[i] += (w[i]*v[i+1] - w[i+1]*(2.0*v[i]+v[i+1]))/6.0;
506 ahwv[i] += (w[i-1]*(v[i-1]+2.0*v[i]) - w[i]*v[i-1])/6.0;
512 const std::vector<Real> &w,
513 const std::vector<Real> &v,
514 const std::vector<Real> &u,
515 const std::vector<Real> &z) {
516 ahwv.assign(u.size(),0.0);
519 const std::vector<Real> &w,
520 const std::vector<Real> &v,
521 const std::vector<Real> &u,
522 const std::vector<Real> &z) {
523 ahwv.assign(z.size(),0.0);
526 const std::vector<Real> &w,
527 const std::vector<Real> &v,
528 const std::vector<Real> &u,
529 const std::vector<Real> &z) {
530 ahwv.assign(z.size(),0.0);
537 Teuchos::RCP<std::vector<Real> >
vec_;
538 Teuchos::RCP<BurgersFEM<Real> >
fem_;
545 : vec_(vec), fem_(fem), dual_vec_(Teuchos::null) {}
549 const std::vector<Real>& xval = *ex.
getVector();
550 std::copy(xval.begin(),xval.end(),vec_->begin());
555 const std::vector<Real>& xval = *ex.
getVector();
556 unsigned dimension = vec_->size();
557 for (
unsigned i=0; i<dimension; i++) {
558 (*vec_)[i] += xval[i];
563 unsigned dimension = vec_->size();
564 for (
unsigned i=0; i<dimension; i++) {
571 const std::vector<Real>& xval = *ex.
getVector();
572 return fem_->compute_L2_dot(xval,*vec_);
577 val = std::sqrt( dot(*
this) );
581 Teuchos::RCP<ROL::Vector<Real> >
clone()
const {
582 return Teuchos::rcp(
new L2VectorPrimal( Teuchos::rcp(
new std::vector<Real>(vec_->size(),0.0)),fem_));
585 Teuchos::RCP<const std::vector<Real> >
getVector()
const {
593 Teuchos::RCP<ROL::Vector<Real> >
basis(
const int i )
const {
594 Teuchos::RCP<L2VectorPrimal> e
595 = Teuchos::rcp(
new L2VectorPrimal( Teuchos::rcp(
new std::vector<Real>(vec_->size(),0.0)),fem_));
596 (*e->getVector())[i] = 1.0;
606 Teuchos::rcp(
new std::vector<Real>(*vec_)),fem_));
608 fem_->apply_mass(*(Teuchos::rcp_const_cast<std::vector<Real> >(dual_vec_->getVector())),*vec_);
617 Teuchos::RCP<std::vector<Real> >
vec_;
618 Teuchos::RCP<BurgersFEM<Real> >
fem_;
625 : vec_(vec), fem_(fem), dual_vec_(Teuchos::null) {}
629 const std::vector<Real>& xval = *ex.
getVector();
630 std::copy(xval.begin(),xval.end(),vec_->begin());
635 const std::vector<Real>& xval = *ex.
getVector();
636 unsigned dimension = vec_->size();
637 for (
unsigned i=0; i<dimension; i++) {
638 (*vec_)[i] += xval[i];
643 unsigned dimension = vec_->size();
644 for (
unsigned i=0; i<dimension; i++) {
651 const std::vector<Real>& xval = *ex.
getVector();
652 unsigned dimension = vec_->size();
653 std::vector<Real> Mx(dimension,0.0);
654 fem_->apply_inverse_mass(Mx,xval);
656 for (
unsigned i = 0; i < dimension; i++) {
657 val += Mx[i]*(*vec_)[i];
664 val = std::sqrt( dot(*
this) );
668 Teuchos::RCP<ROL::Vector<Real> >
clone()
const {
669 return Teuchos::rcp(
new L2VectorDual( Teuchos::rcp(
new std::vector<Real>(vec_->size(),0.0)),fem_));
672 Teuchos::RCP<const std::vector<Real> >
getVector()
const {
680 Teuchos::RCP<ROL::Vector<Real> >
basis(
const int i )
const {
681 Teuchos::RCP<L2VectorDual> e
682 = Teuchos::rcp(
new L2VectorDual( Teuchos::rcp(
new std::vector<Real>(vec_->size(),0.0)),fem_));
683 (*e->getVector())[i] = 1.0;
693 Teuchos::rcp(
new std::vector<Real>(*vec_)),fem_));
695 fem_->apply_inverse_mass(*(Teuchos::rcp_const_cast<std::vector<Real> >(dual_vec_->getVector())),*vec_);
704 Teuchos::RCP<std::vector<Real> >
vec_;
705 Teuchos::RCP<BurgersFEM<Real> >
fem_;
712 : vec_(vec), fem_(fem), dual_vec_(Teuchos::null) {}
716 const std::vector<Real>& xval = *ex.
getVector();
717 std::copy(xval.begin(),xval.end(),vec_->begin());
722 const std::vector<Real>& xval = *ex.
getVector();
723 unsigned dimension = vec_->size();
724 for (
unsigned i=0; i<dimension; i++) {
725 (*vec_)[i] += xval[i];
730 unsigned dimension = vec_->size();
731 for (
unsigned i=0; i<dimension; i++) {
738 const std::vector<Real>& xval = *ex.
getVector();
739 return fem_->compute_H1_dot(xval,*vec_);
744 val = std::sqrt( dot(*
this) );
748 Teuchos::RCP<ROL::Vector<Real> >
clone()
const {
749 return Teuchos::rcp(
new H1VectorPrimal( Teuchos::rcp(
new std::vector<Real>(vec_->size(),0.0)),fem_));
752 Teuchos::RCP<const std::vector<Real> >
getVector()
const {
760 Teuchos::RCP<ROL::Vector<Real> >
basis(
const int i )
const {
761 Teuchos::RCP<H1VectorPrimal> e
762 = Teuchos::rcp(
new H1VectorPrimal( Teuchos::rcp(
new std::vector<Real>(vec_->size(),0.0)),fem_));
763 (*e->getVector())[i] = 1.0;
773 Teuchos::rcp(
new std::vector<Real>(*vec_)),fem_));
775 fem_->apply_H1(*(Teuchos::rcp_const_cast<std::vector<Real> >(dual_vec_->getVector())),*vec_);
784 Teuchos::RCP<std::vector<Real> >
vec_;
785 Teuchos::RCP<BurgersFEM<Real> >
fem_;
792 : vec_(vec), fem_(fem), dual_vec_(Teuchos::null) {}
796 const std::vector<Real>& xval = *ex.
getVector();
797 std::copy(xval.begin(),xval.end(),vec_->begin());
802 const std::vector<Real>& xval = *ex.
getVector();
803 unsigned dimension = vec_->size();
804 for (
unsigned i=0; i<dimension; i++) {
805 (*vec_)[i] += xval[i];
810 unsigned dimension = vec_->size();
811 for (
unsigned i=0; i<dimension; i++) {
818 const std::vector<Real>& xval = *ex.
getVector();
819 unsigned dimension = vec_->size();
820 std::vector<Real> Mx(dimension,0.0);
821 fem_->apply_inverse_H1(Mx,xval);
823 for (
unsigned i = 0; i < dimension; i++) {
824 val += Mx[i]*(*vec_)[i];
831 val = std::sqrt( dot(*
this) );
835 Teuchos::RCP<ROL::Vector<Real> >
clone()
const {
836 return Teuchos::rcp(
new H1VectorDual( Teuchos::rcp(
new std::vector<Real>(vec_->size(),0.0)),fem_));
839 Teuchos::RCP<const std::vector<Real> >
getVector()
const {
847 Teuchos::RCP<ROL::Vector<Real> >
basis(
const int i )
const {
848 Teuchos::RCP<H1VectorDual> e
849 = Teuchos::rcp(
new H1VectorDual( Teuchos::rcp(
new std::vector<Real>(vec_->size(),0.0)),fem_));
850 (*e->getVector())[i] = 1.0;
860 Teuchos::rcp(
new std::vector<Real>(*vec_)),fem_));
862 fem_->apply_inverse_H1(*(Teuchos::rcp_const_cast<std::vector<Real> >(dual_vec_->getVector())),*vec_);
881 Teuchos::RCP<BurgersFEM<Real> >
fem_;
886 const bool useHessian =
true)
887 : fem_(fem), useHessian_(useHessian) {}
891 Teuchos::RCP<std::vector<Real> > cp =
892 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<PrimalConstraintVector>(c)).getVector());
893 Teuchos::RCP<const std::vector<Real> > up =
894 (Teuchos::dyn_cast<PrimalStateVector>(
const_cast<ROL::Vector<Real> &
>(u))).getVector();
895 Teuchos::RCP<const std::vector<Real> > zp =
896 (Teuchos::dyn_cast<PrimalControlVector>(
const_cast<ROL::Vector<Real> &
>(z))).getVector();
898 fem_->compute_residual(*cp,*up,*zp);
903 Teuchos::RCP<std::vector<Real> > jvp =
904 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<PrimalConstraintVector>(jv)).getVector());
905 Teuchos::RCP<const std::vector<Real> > vp =
906 (Teuchos::dyn_cast<PrimalStateVector>(
const_cast<ROL::Vector<Real> &
>(v))).getVector();
907 Teuchos::RCP<const std::vector<Real> > up =
908 (Teuchos::dyn_cast<PrimalStateVector>(
const_cast<ROL::Vector<Real> &
>(u))).getVector();
909 Teuchos::RCP<const std::vector<Real> > zp =
910 (Teuchos::dyn_cast<PrimalControlVector>(
const_cast<ROL::Vector<Real> &
>(z))).getVector();
912 fem_->apply_pde_jacobian(*jvp,*vp,*up,*zp);
917 Teuchos::RCP<std::vector<Real> > jvp =
918 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<PrimalConstraintVector>(jv)).getVector());
919 Teuchos::RCP<const std::vector<Real> > vp =
920 (Teuchos::dyn_cast<PrimalControlVector>(
const_cast<ROL::Vector<Real> &
>(v))).getVector();
921 Teuchos::RCP<const std::vector<Real> > up =
922 (Teuchos::dyn_cast<PrimalStateVector>(
const_cast<ROL::Vector<Real> &
>(u))).getVector();
923 Teuchos::RCP<const std::vector<Real> > zp =
924 (Teuchos::dyn_cast<PrimalControlVector>(
const_cast<ROL::Vector<Real> &
>(z))).getVector();
926 fem_->apply_control_jacobian(*jvp,*vp,*up,*zp);
931 Teuchos::RCP<std::vector<Real> > ijvp =
932 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<PrimalStateVector>(ijv)).getVector());
933 Teuchos::RCP<const std::vector<Real> > vp =
934 (Teuchos::dyn_cast<PrimalConstraintVector>(
const_cast<ROL::Vector<Real> &
>(v))).getVector();
935 Teuchos::RCP<const std::vector<Real> > up =
936 (Teuchos::dyn_cast<PrimalStateVector>(
const_cast<ROL::Vector<Real> &
>(u))).getVector();
937 Teuchos::RCP<const std::vector<Real> > zp =
938 (Teuchos::dyn_cast<PrimalControlVector>(
const_cast<ROL::Vector<Real> &
>(z))).getVector();
940 fem_->apply_inverse_pde_jacobian(*ijvp,*vp,*up,*zp);
945 Teuchos::RCP<std::vector<Real> > jvp =
946 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<DualStateVector>(ajv)).getVector());
947 Teuchos::RCP<const std::vector<Real> > vp =
948 (Teuchos::dyn_cast<DualConstraintVector>(
const_cast<ROL::Vector<Real> &
>(v))).getVector();
949 Teuchos::RCP<const std::vector<Real> > up =
950 (Teuchos::dyn_cast<PrimalStateVector>(
const_cast<ROL::Vector<Real> &
>(u))).getVector();
951 Teuchos::RCP<const std::vector<Real> > zp =
952 (Teuchos::dyn_cast<PrimalControlVector>(
const_cast<ROL::Vector<Real> &
>(z))).getVector();
954 fem_->apply_adjoint_pde_jacobian(*jvp,*vp,*up,*zp);
959 Teuchos::RCP<std::vector<Real> > jvp =
960 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<DualControlVector>(jv)).getVector());
961 Teuchos::RCP<const std::vector<Real> > vp =
962 (Teuchos::dyn_cast<DualConstraintVector>(
const_cast<ROL::Vector<Real> &
>(v))).getVector();
963 Teuchos::RCP<const std::vector<Real> > up =
964 (Teuchos::dyn_cast<PrimalStateVector>(
const_cast<ROL::Vector<Real> &
>(u))).getVector();
965 Teuchos::RCP<const std::vector<Real> > zp =
966 (Teuchos::dyn_cast<PrimalControlVector>(
const_cast<ROL::Vector<Real> &
>(z))).getVector();
968 fem_->apply_adjoint_control_jacobian(*jvp,*vp,*up,*zp);
973 Teuchos::RCP<std::vector<Real> > iajvp =
974 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<DualConstraintVector>(iajv)).getVector());
975 Teuchos::RCP<const std::vector<Real> > vp =
976 (Teuchos::dyn_cast<DualStateVector>(
const_cast<ROL::Vector<Real> &
>(v))).getVector();
977 Teuchos::RCP<const std::vector<Real> > up =
978 (Teuchos::dyn_cast<PrimalStateVector>(
const_cast<ROL::Vector<Real> &
>(u))).getVector();
979 Teuchos::RCP<const std::vector<Real> > zp =
980 (Teuchos::dyn_cast<PrimalControlVector>(
const_cast<ROL::Vector<Real> &
>(z))).getVector();
982 fem_->apply_inverse_adjoint_pde_jacobian(*iajvp,*vp,*up,*zp);
988 Teuchos::RCP<std::vector<Real> > ahwvp =
989 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<DualStateVector>(ahwv)).getVector());
990 Teuchos::RCP<const std::vector<Real> > wp =
991 (Teuchos::dyn_cast<DualConstraintVector>(
const_cast<ROL::Vector<Real> &
>(w))).getVector();
992 Teuchos::RCP<const std::vector<Real> > vp =
993 (Teuchos::dyn_cast<PrimalStateVector>(
const_cast<ROL::Vector<Real> &
>(v))).getVector();
994 Teuchos::RCP<const std::vector<Real> > up =
995 (Teuchos::dyn_cast<PrimalStateVector>(
const_cast<ROL::Vector<Real> &
>(u))).getVector();
996 Teuchos::RCP<const std::vector<Real> > zp =
997 (Teuchos::dyn_cast<PrimalControlVector>(
const_cast<ROL::Vector<Real> &
>(z))).getVector();
999 fem_->apply_adjoint_pde_hessian(*ahwvp,*wp,*vp,*up,*zp);
1008 if ( useHessian_ ) {
1009 Teuchos::RCP<std::vector<Real> > ahwvp =
1010 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<DualControlVector>(ahwv)).getVector());
1011 Teuchos::RCP<const std::vector<Real> > wp =
1012 (Teuchos::dyn_cast<DualConstraintVector>(
const_cast<ROL::Vector<Real> &
>(w))).getVector();
1013 Teuchos::RCP<const std::vector<Real> > vp =
1014 (Teuchos::dyn_cast<PrimalStateVector>(
const_cast<ROL::Vector<Real> &
>(v))).getVector();
1015 Teuchos::RCP<const std::vector<Real> > up =
1016 (Teuchos::dyn_cast<PrimalStateVector>(
const_cast<ROL::Vector<Real> &
>(u))).getVector();
1017 Teuchos::RCP<const std::vector<Real> > zp =
1018 (Teuchos::dyn_cast<PrimalControlVector>(
const_cast<ROL::Vector<Real> &
>(z))).getVector();
1020 fem_->apply_adjoint_control_pde_hessian(*ahwvp,*wp,*vp,*up,*zp);
1028 if ( useHessian_ ) {
1029 Teuchos::RCP<std::vector<Real> > ahwvp =
1030 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<DualStateVector>(ahwv)).getVector());
1031 Teuchos::RCP<const std::vector<Real> > wp =
1032 (Teuchos::dyn_cast<DualConstraintVector>(
const_cast<ROL::Vector<Real> &
>(w))).getVector();
1033 Teuchos::RCP<const std::vector<Real> > vp =
1034 (Teuchos::dyn_cast<PrimalControlVector>(
const_cast<ROL::Vector<Real> &
>(v))).getVector();
1035 Teuchos::RCP<const std::vector<Real> > up =
1036 (Teuchos::dyn_cast<PrimalStateVector>(
const_cast<ROL::Vector<Real> &
>(u))).getVector();
1037 Teuchos::RCP<const std::vector<Real> > zp =
1038 (Teuchos::dyn_cast<PrimalControlVector>(
const_cast<ROL::Vector<Real> &
>(z))).getVector();
1040 fem_->apply_adjoint_pde_control_hessian(*ahwvp,*wp,*vp,*up,*zp);
1048 if ( useHessian_ ) {
1049 Teuchos::RCP<std::vector<Real> > ahwvp =
1050 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<DualControlVector>(ahwv)).getVector());
1051 Teuchos::RCP<const std::vector<Real> > wp =
1052 (Teuchos::dyn_cast<DualConstraintVector>(
const_cast<ROL::Vector<Real> &
>(w))).getVector();
1053 Teuchos::RCP<const std::vector<Real> > vp =
1054 (Teuchos::dyn_cast<PrimalControlVector>(
const_cast<ROL::Vector<Real> &
>(v))).getVector();
1055 Teuchos::RCP<const std::vector<Real> > up =
1056 (Teuchos::dyn_cast<PrimalStateVector>(
const_cast<ROL::Vector<Real> &
>(u))).getVector();
1057 Teuchos::RCP<const std::vector<Real> > zp =
1058 (Teuchos::dyn_cast<PrimalControlVector>(
const_cast<ROL::Vector<Real> &
>(z))).getVector();
1060 fem_->apply_adjoint_control_hessian(*ahwvp,*wp,*vp,*up,*zp);
BurgersFEM(int nx=128, Real nl=1.0, Real cH1=1.0, Real cL2=1.0)
L2VectorPrimal< Real > PrimalControlVector
Real compute_H1_dot(const std::vector< Real > &x, const std::vector< Real > &y) const
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 ...
Real norm() const
Returns where .
void apply_mass(std::vector< Real > &Mu, const std::vector< Real > &u) const
Teuchos::RCP< const std::vector< Real > > getVector() const
void axpy(std::vector< Real > &out, const Real a, const std::vector< Real > &x, const std::vector< Real > &y) const
Teuchos::RCP< std::vector< Real > > getVector()
EqualityConstraint_BurgersControl(const Teuchos::RCP< BurgersFEM< Real > > &fem, const bool useHessian=true)
Teuchos::RCP< ROL::Vector< Real > > basis(const int i) const
Return i-th basis vector.
Real dot(const ROL::Vector< Real > &x) const
Compute where .
void apply_adjoint_pde_control_hessian(std::vector< Real > &ahwv, const std::vector< Real > &w, const std::vector< Real > &v, const std::vector< Real > &u, const std::vector< Real > &z)
Teuchos::RCP< const std::vector< Real > > getVector() const
int dimension() const
Return dimension of the vector space.
Teuchos::RCP< std::vector< Real > > vec_
void plus(const ROL::Vector< Real > &x)
Compute , where .
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 .
void apply_adjoint_pde_hessian(std::vector< Real > &ahwv, const std::vector< Real > &w, const std::vector< Real > &v, const std::vector< Real > &u, const std::vector< Real > &z) const
L2VectorDual(const Teuchos::RCP< std::vector< Real > > &vec, const Teuchos::RCP< BurgersFEM< Real > > &fem)
Teuchos::RCP< ROL::Vector< Real > > basis(const int i) const
Return i-th basis vector.
Teuchos::RCP< H1VectorPrimal< Real > > dual_vec_
Real norm() const
Returns where .
Contains definitions of custom data types in ROL.
Real dot(const ROL::Vector< Real > &x) const
Compute where .
void plus(const ROL::Vector< Real > &x)
Compute , where .
Teuchos::RCP< const std::vector< Real > > getVector() const
Real compute_H1_norm(const std::vector< Real > &r) const
Teuchos::RCP< std::vector< Real > > vec_
Real norm() const
Returns where .
Teuchos::RCP< std::vector< Real > > getVector()
Teuchos::RCP< BurgersFEM< Real > > fem_
Teuchos::RCP< std::vector< Real > > getVector()
virtual void zero()
Set to zero vector.
Defines the linear algebra or vector space interface.
Teuchos::RCP< H1VectorDual< Real > > dual_vec_
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 adjoint of the partial constraint Hessian at , , to vector in direction ...
void scale(const Real alpha)
Compute where .
Defines the equality constraint operator interface for simulation-based optimization.
L2VectorDual< Real > DualControlVector
Real dot(const ROL::Vector< Real > &x) const
Compute where .
void plus(const ROL::Vector< Real > &x)
Compute , where .
Real mesh_spacing(void) const
void applyInverseAdjointJacobian_1(ROL::Vector< Real > &iajv, 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 ...
int dimension() const
Return dimension of the vector space.
void test_inverse_mass(std::ostream &outStream=std::cout)
void apply_pde_jacobian(std::vector< Real > &jv, const std::vector< Real > &v, const std::vector< Real > &u, const std::vector< Real > &z) const
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...
Teuchos::RCP< BurgersFEM< Real > > fem_
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 .
Teuchos::RCP< ROL::Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
void scale(const Real alpha)
Compute where .
void compute_pde_jacobian(std::vector< Real > &dl, std::vector< Real > &d, std::vector< Real > &du, const std::vector< Real > &u) const
H1VectorPrimal< Real > PrimalStateVector
Teuchos::RCP< ROL::Vector< Real > > basis(const int i) const
Return i-th basis vector.
Real compute_L2_dot(const std::vector< Real > &x, const std::vector< Real > &y) const
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) const
Teuchos::RCP< L2VectorPrimal< Real > > dual_vec_
Teuchos::RCP< BurgersFEM< Real > > fem_
void apply_inverse_mass(std::vector< Real > &Mu, const std::vector< Real > &u) const
H1VectorDual< Real > PrimalConstraintVector
void update(std::vector< Real > &u, const std::vector< Real > &s, const Real alpha=1.0) const
void scale(const Real alpha)
Compute where .
H1VectorPrimal< Real > DualConstraintVector
Teuchos::RCP< std::vector< Real > > vec_
void apply_adjoint_control_hessian(std::vector< Real > &ahwv, const std::vector< Real > &w, const std::vector< Real > &v, const std::vector< Real > &u, const std::vector< Real > &z)
void compute_residual(std::vector< Real > &r, const std::vector< Real > &u, const std::vector< Real > &z) const
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 .
H1VectorDual< Real > DualStateVector
void set_problem_data(const Real nu, const Real f, const Real u0, const Real u1)
void apply_control_jacobian(std::vector< Real > &jv, const std::vector< Real > &v, const std::vector< Real > &u, const std::vector< Real > &z) const
void test_inverse_H1(std::ostream &outStream=std::cout)
void value(ROL::Vector< Real > &c, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Evaluate the constraint operator at .
const ROL::Vector< Real > & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
H1VectorDual(const Teuchos::RCP< std::vector< Real > > &vec, const Teuchos::RCP< BurgersFEM< Real > > &fem)
void apply_adjoint_control_jacobian(std::vector< Real > &jv, const std::vector< Real > &v, const std::vector< Real > &u, const std::vector< Real > &z) const
Real dot(const ROL::Vector< Real > &x) const
Compute where .
Teuchos::RCP< BurgersFEM< Real > > fem_
Teuchos::RCP< L2VectorDual< Real > > dual_vec_
L2VectorPrimal(const Teuchos::RCP< std::vector< Real > > &vec, const Teuchos::RCP< BurgersFEM< Real > > &fem)
void apply_inverse_pde_jacobian(std::vector< Real > &ijv, const std::vector< Real > &v, const std::vector< Real > &u, const std::vector< Real > &z) const
Teuchos::RCP< ROL::Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
Teuchos::RCP< std::vector< Real > > vec_
void apply_H1(std::vector< Real > &Mu, const std::vector< Real > &u) const
void plus(const ROL::Vector< Real > &x)
Compute , where .
void apply_inverse_H1(std::vector< Real > &Mu, const std::vector< Real > &u) const
void apply_adjoint_pde_jacobian(std::vector< Real > &ajv, const std::vector< Real > &v, const std::vector< Real > &u, const std::vector< Real > &z) const
void apply_inverse_adjoint_pde_jacobian(std::vector< Real > &iajv, const std::vector< Real > &v, const std::vector< Real > &u, const std::vector< Real > &z) const
void apply_adjoint_control_pde_hessian(std::vector< Real > &ahwv, const std::vector< Real > &w, const std::vector< Real > &v, const std::vector< Real > &u, const std::vector< Real > &z)
Teuchos::RCP< BurgersFEM< Real > > fem_
int dimension() const
Return dimension of the vector space.
const ROL::Vector< Real > & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
Teuchos::RCP< ROL::Vector< Real > > basis(const int i) const
Return i-th basis vector.
H1VectorPrimal(const Teuchos::RCP< std::vector< Real > > &vec, const Teuchos::RCP< BurgersFEM< Real > > &fem)
Real norm() const
Returns where .
void scale(std::vector< Real > &u, const Real alpha=0.0) const
Teuchos::RCP< std::vector< Real > > getVector()
const ROL::Vector< Real > & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
Teuchos::RCP< ROL::Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
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 ...
Teuchos::RCP< const std::vector< Real > > getVector() const
void scale(const Real alpha)
Compute where .
int dimension() const
Return dimension of the vector space.
const ROL::Vector< Real > & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
Real compute_L2_norm(const std::vector< Real > &r) const
Teuchos::RCP< ROL::Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.