65 return std::sqrt(this->
dot(r,r));
68 Real
dot(
const std::vector<Real> &x,
const std::vector<Real> &y) {
70 Real c = (((int)x.size()==this->
nx_) ? 4.0 : 2.0);
71 for (
unsigned i=0; i<x.size(); i++) {
73 ip += this->
dx_/6.0*(c*x[i] + x[i+1])*y[i];
75 else if ( i == x.size()-1 ) {
76 ip += this->
dx_/6.0*(x[i-1] + c*x[i])*y[i];
79 ip += this->
dx_/6.0*(x[i-1] + 4.0*x[i] + x[i+1])*y[i];
87 void update(std::vector<Real> &u,
const std::vector<Real> &s,
const Real alpha=1.0) {
88 for (
unsigned i=0; i<u.size(); i++) {
93 void scale(std::vector<Real> &u,
const Real alpha=0.0) {
94 for (
unsigned i=0; i<u.size(); i++) {
100 const std::vector<Real> &z) {
102 r.resize(this->
nx_,0.0);
103 for (
int i=0; i<this->
nx_; i++) {
106 r[i] = this->
nu_/this->
dx_*(2.0*u[i]-u[i+1]);
108 else if (i==this->nx_-1) {
109 r[i] = this->
nu_/this->
dx_*(2.0*u[i]-u[i-1]);
112 r[i] = this->
nu_/this->
dx_*(2.0*u[i]-u[i-1]-u[i+1]);
116 r[i] += u[i+1]*(u[i]+u[i+1])/6.0;
119 r[i] -= u[i-1]*(u[i-1]+u[i])/6.0;
122 r[i] -= this->
dx_/6.0*(z[i]+4.0*z[i+1]+z[i+2]);
124 r[i] -= this->
dx_*this->
f_;
127 r[0] -= this->
u0_*u[ 0]/6.0 + this->
u0_*this->
u0_/6.0 + this->
nu_*this->
u0_/this->
dx_;
128 r[this->nx_-1] += this->
u1_*u[this->nx_-1]/6.0 + this->
u1_*this->
u1_/6.0 - this->
nu_*this->
u1_/this->
dx_;
132 const std::vector<Real> &u) {
135 d.resize(this->
nx_,this->
nu_*2.0/this->
dx_);
137 dl.resize(this->
nx_-1,-this->
nu_/this->
dx_);
139 du.resize(this->
nx_-1,-this->
nu_/this->
dx_);
141 for (
int i=0; i<this->
nx_; i++) {
143 dl[i] += (-2.0*u[i]-u[i+1])/6.0;
148 du[i-1] += (u[i-1]+2.0*u[i])/6.0;
152 d[0] -= this->
u0_/6.0;
153 d[this->nx_-1] += this->
u1_/6.0;
156 void linear_solve(std::vector<Real> &u, std::vector<Real> &dl, std::vector<Real> &d, std::vector<Real> &du,
157 const std::vector<Real> &r,
const bool transpose =
false) {
158 u.assign(r.begin(),r.end());
160 Teuchos::LAPACK<int,Real> lp;
161 std::vector<Real> du2(this->
nx_-2,0.0);
162 std::vector<int> ipiv(this->
nx_,0);
166 lp.GTTRF(this->
nx_,&dl[0],&d[0],&du[0],&du2[0],&ipiv[0],&info);
171 lp.GTTRS(trans,this->
nx_,nhrs,&dl[0],&d[0],&du[0],&du2[0],&ipiv[0],&u[0],ldb,&info);
178 dx_ = 1.0/((Real)nx+1.0);
183 ROL::Ptr<std::vector<Real> > cp =
185 ROL::Ptr<const std::vector<Real> > up =
187 ROL::Ptr<const std::vector<Real> > zp =
194 ROL::Ptr<std::vector<Real> > jvp =
196 ROL::Ptr<const std::vector<Real> > vp =
198 ROL::Ptr<const std::vector<Real> > up =
200 ROL::Ptr<const std::vector<Real> > zp =
203 for (
int i = 0; i < this->
nx_; i++) {
204 (*jvp)[i] = this->
nu_/this->
dx_*2.0*(*vp)[i];
206 (*jvp)[i] += -this->
nu_/this->
dx_*(*vp)[i-1]
207 -(*up)[i-1]/6.0*(*vp)[i]
208 -((*up)[i]+2.0*(*up)[i-1])/6.0*(*vp)[i-1];
210 if ( i < this->nx_-1 ) {
211 (*jvp)[i] += -this->
nu_/this->
dx_*(*vp)[i+1]
212 +(*up)[i+1]/6.0*(*vp)[i]
213 +((*up)[i]+2.0*(*up)[i+1])/6.0*(*vp)[i+1];
216 (*jvp)[0] -= this->
u0_/6.0*(*vp)[0];
217 (*jvp)[this->nx_-1] += this->
u1_/6.0*(*vp)[this->nx_-1];
222 ROL::Ptr<std::vector<Real> > jvp =
224 ROL::Ptr<const std::vector<Real> > vp =
226 ROL::Ptr<const std::vector<Real> > up =
228 ROL::Ptr<const std::vector<Real> > zp =
230 for (
int i=0; i<this->
nx_; i++) {
232 (*jvp)[i] = -this->
dx_/6.0*((*vp)[i]+4.0*(*vp)[i+1]+(*vp)[i+2]);
238 ROL::Ptr<std::vector<Real> > ijvp =
240 ROL::Ptr<const std::vector<Real> > vp =
242 ROL::Ptr<const std::vector<Real> > up =
244 ROL::Ptr<const std::vector<Real> > zp =
247 std::vector<Real> d(this->
nx_,0.0);
248 std::vector<Real> dl(this->
nx_-1,0.0);
249 std::vector<Real> du(this->
nx_-1,0.0);
257 ROL::Ptr<std::vector<Real> > jvp =
259 ROL::Ptr<const std::vector<Real> > vp =
261 ROL::Ptr<const std::vector<Real> > up =
263 ROL::Ptr<const std::vector<Real> > zp =
266 for (
int i = 0; i < this->
nx_; i++) {
267 (*jvp)[i] = this->
nu_/this->
dx_*2.0*(*vp)[i];
269 (*jvp)[i] += -this->
nu_/this->
dx_*(*vp)[i-1]
270 -(*up)[i-1]/6.0*(*vp)[i]
271 +((*up)[i-1]+2.0*(*up)[i])/6.0*(*vp)[i-1];
273 if ( i < this->nx_-1 ) {
274 (*jvp)[i] += -this->
nu_/this->
dx_*(*vp)[i+1]
275 +(*up)[i+1]/6.0*(*vp)[i]
276 -((*up)[i+1]+2.0*(*up)[i])/6.0*(*vp)[i+1];
279 (*jvp)[0] -= this->
u0_/6.0*(*vp)[0];
280 (*jvp)[this->nx_-1] += this->
u1_/6.0*(*vp)[this->nx_-1];
285 ROL::Ptr<std::vector<Real> > jvp =
287 ROL::Ptr<const std::vector<Real> > vp =
289 ROL::Ptr<const std::vector<Real> > up =
291 ROL::Ptr<const std::vector<Real> > zp =
293 for (
int i=0; i<this->
nx_+2; i++) {
295 (*jvp)[i] = -this->
dx_/6.0*(*vp)[i];
298 (*jvp)[i] = -this->
dx_/6.0*(4.0*(*vp)[i-1]+(*vp)[i]);
300 else if ( i == this->nx_ ) {
301 (*jvp)[i] = -this->
dx_/6.0*(4.0*(*vp)[i-1]+(*vp)[i-2]);
303 else if ( i == this->nx_+1 ) {
304 (*jvp)[i] = -this->
dx_/6.0*(*vp)[i-2];
307 (*jvp)[i] = -this->
dx_/6.0*((*vp)[i-2]+4.0*(*vp)[i-1]+(*vp)[i]);
314 ROL::Ptr<std::vector<Real> > iajvp =
316 ROL::Ptr<const std::vector<Real> > vp =
318 ROL::Ptr<const std::vector<Real> > up =
321 std::vector<Real> d(this->
nx_,0.0);
322 std::vector<Real> du(this->
nx_-1,0.0);
323 std::vector<Real> dl(this->
nx_-1,0.0);
331 ROL::Ptr<std::vector<Real> > ahwvp =
333 ROL::Ptr<const std::vector<Real> > wp =
335 ROL::Ptr<const std::vector<Real> > vp =
337 ROL::Ptr<const std::vector<Real> > up =
339 ROL::Ptr<const std::vector<Real> > zp =
341 for (
int i=0; i<this->
nx_; i++) {
345 (*ahwvp)[i] += ((*wp)[i]*(*vp)[i+1] - (*wp)[i+1]*(2.0*(*vp)[i]+(*vp)[i+1]))/6.0;
348 (*ahwvp)[i] += ((*wp)[i-1]*((*vp)[i-1]+2.0*(*vp)[i]) - (*wp)[i]*(*vp)[i-1])/6.0;
385 case 1: val = ((x<0.5) ? 1.0 : 0.0);
break;
386 case 2: val = 1.0;
break;
387 case 3: val = std::abs(std::sin(8.0*M_PI*x));
break;
388 case 4: val = std::exp(-0.5*(x-0.5)*(x-0.5));
break;
393 Real
dot(
const std::vector<Real> &x,
const std::vector<Real> &y) {
395 Real c = (((int)x.size()==this->
nx_) ? 4.0 : 2.0);
396 for (
unsigned i=0; i<x.size(); i++) {
398 ip += this->
dx_/6.0*(c*x[i] + x[i+1])*y[i];
400 else if ( i == x.size()-1 ) {
401 ip += this->
dx_/6.0*(x[i-1] + c*x[i])*y[i];
404 ip += this->
dx_/6.0*(x[i-1] + 4.0*x[i] + x[i+1])*y[i];
410 void apply_mass(std::vector<Real> &Mu,
const std::vector<Real> &u ) {
411 Mu.resize(u.size(),0.0);
412 Real c = (((int)u.size()==this->
nx_) ? 4.0 : 2.0);
413 for (
unsigned i=0; i<u.size(); i++) {
415 Mu[i] = this->
dx_/6.0*(c*u[i] + u[i+1]);
417 else if ( i == u.size()-1 ) {
418 Mu[i] = this->
dx_/6.0*(u[i-1] + c*u[i]);
421 Mu[i] = this->
dx_/6.0*(u[i-1] + 4.0*u[i] + u[i+1]);
432 dx_ = 1.0/((Real)nx+1.0);
436 ROL::Ptr<const std::vector<Real> > up =
438 ROL::Ptr<const std::vector<Real> > zp =
441 Real res1 = 0.0, res2 = 0.0, res3 = 0.0;
442 Real valu = 0.0, valz = this->
dot(*zp,*zp);
443 for (
int i=0; i<this->
nx_; i++) {
447 valu += this->
dx_/6.0*(4.0*res1 + res2)*res1;
449 else if ( i == this->nx_-1 ) {
452 valu += this->dx_/6.0*(res1 + 4.0*res2)*res2;
458 valu += this->dx_/6.0*(res1 + 4.0*res2 + res3)*res2;
461 return 0.5*(valu + this->
alpha_*valz);
466 ROL::Ptr<std::vector<Real> > gup = ROL::constPtrCast<std::vector<Real> >(
469 ROL::Ptr<const std::vector<Real> > up =
471 ROL::Ptr<const std::vector<Real> > zp =
474 std::vector<Real> diff(this->
nx_,0.0);
475 for (
int i=0; i<this->
nx_; i++) {
483 ROL::Ptr<std::vector<Real> > gzp = ROL::constPtrCast<std::vector<Real> >(
486 ROL::Ptr<const std::vector<Real> > up =
488 ROL::Ptr<const std::vector<Real> > zp =
491 for (
int i=0; i<this->
nx_+2; i++) {
493 (*gzp)[i] = this->
alpha_*this->
dx_/6.0*(2.0*(*zp)[i]+(*zp)[i+1]);
495 else if (i==this->nx_+1) {
496 (*gzp)[i] = this->
alpha_*this->
dx_/6.0*(2.0*(*zp)[i]+(*zp)[i-1]);
499 (*gzp)[i] = this->
alpha_*this->
dx_/6.0*((*zp)[i-1]+4.0*(*zp)[i]+(*zp)[i+1]);
506 ROL::Ptr<std::vector<Real> > hvup =
509 ROL::Ptr<const std::vector<Real> > vup =
527 ROL::Ptr<std::vector<Real> > hvzp =
530 ROL::Ptr<const std::vector<Real> > vzp =
533 for (
int i=0; i<this->
nx_+2; i++) {
535 (*hvzp)[i] = this->
alpha_*this->
dx_/6.0*(2.0*(*vzp)[i]+(*vzp)[i+1]);
537 else if (i==this->nx_+1) {
538 (*hvzp)[i] = this->
alpha_*this->
dx_/6.0*(2.0*(*vzp)[i]+(*vzp)[i-1]);
541 (*hvzp)[i] = this->
alpha_*this->
dx_/6.0*((*vzp)[i-1]+4.0*(*vzp)[i]+(*vzp)[i+1]);
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 > &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 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 ...
virtual void zero()
Set to zero vector.
Defines the linear algebra or vector space interface.
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 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)
Objective_BurgersControl(Real alpha=1.e-4, int nx=128)
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 hessVec_12(ROL::Vector< Real > &hv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Constraint_BurgersControl(int nx=128, Real nu=1.e-2, Real u0=1.0, Real u1=0.0, Real f=0.0)
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)
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 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 ...
Real compute_norm(const std::vector< Real > &r)
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)