44 #include "Teuchos_ParameterList.hpp" 45 #include "Teuchos_XMLParameterListHelpers.hpp" 46 #include "Teuchos_oblackholestream.hpp" 47 #include "Teuchos_LAPACK.hpp" 48 #include "Teuchos_GlobalMPISession.hpp" 49 #include "Teuchos_Comm.hpp" 50 #include "Teuchos_DefaultComm.hpp" 51 #include "Teuchos_CommHelpers.hpp" 75 return std::sqrt(
dot(r,r));
78 Real
dot(
const std::vector<Real> &x,
const std::vector<Real> &y) {
80 Real c = (((int)x.size()==
nx_) ? 4.0 : 2.0);
81 for (
unsigned i=0; i<x.size(); i++) {
83 ip += dx_/6.0*(c*x[i] + x[i+1])*y[i];
85 else if ( i == x.size()-1 ) {
86 ip += dx_/6.0*(x[i-1] + c*x[i])*y[i];
89 ip += dx_/6.0*(x[i-1] + 4.0*x[i] + x[i+1])*y[i];
95 void update(std::vector<Real> &u,
const std::vector<Real> &s,
const Real alpha=1.0) {
96 for (
unsigned i=0; i<u.size(); i++) {
101 void scale(std::vector<Real> &u,
const Real alpha=0.0) {
102 for (
unsigned i=0; i<u.size(); i++) {
108 const std::vector<Real> &z) {
109 r.clear(); r.resize(nx_,0.0);
110 const std::vector<Real> param =
112 Real nu = std::pow(10.0,param[0]-2.0);
113 Real f = param[1]/100.0;
114 Real u0 = 1.0+param[2]/1000.0;
115 Real u1 = param[3]/1000.0;
116 for (
int i=0; i<
nx_; i++) {
119 r[i] = nu/dx_*(2.0*u[i]-u[i+1]);
122 r[i] = nu/dx_*(2.0*u[i]-u[i-1]);
125 r[i] = nu/dx_*(2.0*u[i]-u[i-1]-u[i+1]);
129 r[i] += u[i+1]*(u[i]+u[i+1])/6.0;
132 r[i] -= u[i-1]*(u[i-1]+u[i])/6.0;
135 r[i] -= dx_/6.0*(z[i]+4.0*z[i+1]+z[i+2]);
140 r[ 0] -= u0*u[ 0]/6.0 + u0*u0/6.0 + nu*u0/
dx_;
141 r[nx_-1] += u1*u[nx_-1]/6.0 + u1*u1/6.0 - nu*u1/
dx_;
145 const std::vector<Real> &u) {
146 const std::vector<Real> param =
148 Real nu = std::pow(10.0,param[0]-2.0);
149 Real u0 = 1.0+param[2]/1000.0;
150 Real u1 = param[3]/1000.0;
152 d.clear(); d.resize(nx_,nu*2.0/dx_);
153 dl.clear(); dl.resize(nx_-1,-nu/dx_);
154 du.clear(); du.resize(nx_-1,-nu/dx_);
156 for (
int i=0; i<
nx_; i++) {
158 dl[i] += (-2.0*u[i]-u[i+1])/6.0;
163 du[i-1] += (u[i-1]+2.0*u[i])/6.0;
171 void linear_solve(std::vector<Real> &u, std::vector<Real> &dl, std::vector<Real> &d, std::vector<Real> &du,
172 const std::vector<Real> &r,
const bool transpose =
false) {
173 u.assign(r.begin(),r.end());
175 Teuchos::LAPACK<int,Real> lp;
176 std::vector<Real> du2(nx_-2,0.0);
177 std::vector<int> ipiv(nx_,0);
181 lp.GTTRF(nx_,&dl[0],&d[0],&du[0],&du2[0],&ipiv[0],&info);
186 lp.GTTRS(trans,nx_,nhrs,&dl[0],&d[0],&du[0],&du2[0],&ipiv[0],&u[0],ldb,&info);
195 Teuchos::RCP<std::vector<Real> > cp =
196 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<
ROL::StdVector<Real> >(c)).getVector());
197 Teuchos::RCP<const std::vector<Real> > up =
199 Teuchos::RCP<const std::vector<Real> > zp =
205 Teuchos::RCP<std::vector<Real> > up =
206 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<
ROL::StdVector<Real> >(u)).getVector());
208 up->assign(up->size(),z.
norm()/up->size());
215 Teuchos::RCP<std::vector<Real> > jvp =
216 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<
ROL::StdVector<Real> >(jv)).getVector());
217 Teuchos::RCP<const std::vector<Real> > vp =
219 Teuchos::RCP<const std::vector<Real> > up =
221 Teuchos::RCP<const std::vector<Real> > zp =
223 const std::vector<Real> param =
225 Real nu = std::pow(10.0,param[0]-2.0);
226 Real u0 = 1.0+param[2]/1000.0;
227 Real u1 = param[3]/1000.0;
229 for (
int i = 0; i <
nx_; i++) {
230 (*jvp)[i] = nu/dx_*2.0*(*vp)[i];
232 (*jvp)[i] += -nu/dx_*(*vp)[i-1]
233 -(*up)[i-1]/6.0*(*vp)[i]
234 -((*up)[i]+2.0*(*up)[i-1])/6.0*(*vp)[i-1];
237 (*jvp)[i] += -nu/dx_*(*vp)[i+1]
238 +(*up)[i+1]/6.0*(*vp)[i]
239 +((*up)[i]+2.0*(*up)[i+1])/6.0*(*vp)[i+1];
242 (*jvp)[ 0] -= u0/6.0*(*vp)[0];
243 (*jvp)[nx_-1] += u1/6.0*(*vp)[nx_-1];
248 Teuchos::RCP<std::vector<Real> > jvp =
249 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<
ROL::StdVector<Real> >(jv)).getVector());
250 Teuchos::RCP<const std::vector<Real> > vp =
252 Teuchos::RCP<const std::vector<Real> > up =
254 Teuchos::RCP<const std::vector<Real> > zp =
256 for (
int i=0; i<
nx_; i++) {
258 (*jvp)[i] = -dx_/6.0*((*vp)[i]+4.0*(*vp)[i+1]+(*vp)[i+2]);
264 Teuchos::RCP<std::vector<Real> > ijvp =
265 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<
ROL::StdVector<Real> >(ijv)).getVector());
266 Teuchos::RCP<const std::vector<Real> > vp =
268 Teuchos::RCP<const std::vector<Real> > up =
270 Teuchos::RCP<const std::vector<Real> > zp =
273 std::vector<Real> d(nx_,0.0);
274 std::vector<Real> dl(nx_-1,0.0);
275 std::vector<Real> du(nx_-1,0.0);
283 Teuchos::RCP<std::vector<Real> > jvp =
284 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<
ROL::StdVector<Real> >(ajv)).getVector());
285 Teuchos::RCP<const std::vector<Real> > vp =
287 Teuchos::RCP<const std::vector<Real> > up =
289 Teuchos::RCP<const std::vector<Real> > zp =
291 const std::vector<Real> param =
293 Real nu = std::pow(10.0,param[0]-2.0);
294 Real u0 = 1.0+param[2]/1000.0;
295 Real u1 = param[3]/1000.0;
297 for (
int i = 0; i <
nx_; i++) {
298 (*jvp)[i] = nu/dx_*2.0*(*vp)[i];
300 (*jvp)[i] += -nu/dx_*(*vp)[i-1]
301 -(*up)[i-1]/6.0*(*vp)[i]
302 +((*up)[i-1]+2.0*(*up)[i])/6.0*(*vp)[i-1];
305 (*jvp)[i] += -nu/dx_*(*vp)[i+1]
306 +(*up)[i+1]/6.0*(*vp)[i]
307 -((*up)[i+1]+2.0*(*up)[i])/6.0*(*vp)[i+1];
310 (*jvp)[ 0] -= u0/6.0*(*vp)[0];
311 (*jvp)[nx_-1] += u1/6.0*(*vp)[nx_-1];
316 Teuchos::RCP<std::vector<Real> > jvp =
317 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<
ROL::StdVector<Real> >(jv)).getVector());
318 Teuchos::RCP<const std::vector<Real> > vp =
320 Teuchos::RCP<const std::vector<Real> > up =
322 Teuchos::RCP<const std::vector<Real> > zp =
324 for (
int i=0; i<nx_+2; i++) {
326 (*jvp)[i] = -dx_/6.0*(*vp)[i];
329 (*jvp)[i] = -dx_/6.0*(4.0*(*vp)[i-1]+(*vp)[i]);
331 else if ( i == nx_ ) {
332 (*jvp)[i] = -dx_/6.0*(4.0*(*vp)[i-1]+(*vp)[i-2]);
334 else if ( i == nx_+1 ) {
335 (*jvp)[i] = -dx_/6.0*(*vp)[i-2];
338 (*jvp)[i] = -dx_/6.0*((*vp)[i-2]+4.0*(*vp)[i-1]+(*vp)[i]);
345 Teuchos::RCP<std::vector<Real> > iajvp =
346 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<
ROL::StdVector<Real> >(iajv)).getVector());
347 Teuchos::RCP<const std::vector<Real> > vp =
349 Teuchos::RCP<const std::vector<Real> > up =
352 std::vector<Real> d(nx_,0.0);
353 std::vector<Real> du(nx_-1,0.0);
354 std::vector<Real> dl(nx_-1,0.0);
362 Teuchos::RCP<std::vector<Real> > ahwvp =
363 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<
ROL::StdVector<Real> >(ahwv)).getVector());
364 Teuchos::RCP<const std::vector<Real> > wp =
366 Teuchos::RCP<const std::vector<Real> > vp =
368 Teuchos::RCP<const std::vector<Real> > up =
370 Teuchos::RCP<const std::vector<Real> > zp =
372 for (
int i=0; i<
nx_; i++) {
376 (*ahwvp)[i] += ((*wp)[i]*(*vp)[i+1] - (*wp)[i+1]*(2.0*(*vp)[i]+(*vp)[i+1]))/6.0;
379 (*ahwvp)[i] += ((*wp)[i-1]*((*vp)[i-1]+2.0*(*vp)[i]) - (*wp)[i]*(*vp)[i-1])/6.0;
413 case 1: val = ((x<0.5) ? 1.0 : 0.0);
break;
414 case 2: val = 1.0;
break;
415 case 3: val = std::abs(std::sin(8.0*M_PI*x));
break;
416 case 4: val = std::exp(-0.5*(x-0.5)*(x-0.5));
break;
421 Real
dot(
const std::vector<Real> &x,
const std::vector<Real> &y) {
423 Real c = (((int)x.size()==
nx_) ? 4.0 : 2.0);
424 for (
unsigned i=0; i<x.size(); i++) {
426 ip += dx_/6.0*(c*x[i] + x[i+1])*y[i];
428 else if ( i == x.size()-1 ) {
429 ip += dx_/6.0*(x[i-1] + c*x[i])*y[i];
432 ip += dx_/6.0*(x[i-1] + 4.0*x[i] + x[i+1])*y[i];
438 void apply_mass(std::vector<Real> &Mu,
const std::vector<Real> &u ) {
439 Mu.resize(u.size(),0.0);
440 Real c = (((int)u.size()==
nx_) ? 4.0 : 2.0);
441 for (
unsigned i=0; i<u.size(); i++) {
443 Mu[i] = dx_/6.0*(c*u[i] + u[i+1]);
445 else if ( i == u.size()-1 ) {
446 Mu[i] = dx_/6.0*(u[i-1] + c*u[i]);
449 Mu[i] = dx_/6.0*(u[i-1] + 4.0*u[i] + u[i+1]);
460 dx_ = 1.0/((Real)nx+1.0);
464 Teuchos::RCP<const std::vector<Real> > up =
466 Teuchos::RCP<const std::vector<Real> > zp =
469 Real res1 = 0.0, res2 = 0.0, res3 = 0.0;
470 Real valu = 0.0, valz =
dot(*zp,*zp);
471 for (
int i=0; i<
nx_; i++) {
473 res1 = (*up)[i]-evaluate_target((Real)(i+1)*dx_);
474 res2 = (*up)[i+1]-evaluate_target((Real)(i+2)*dx_);
475 valu += dx_/6.0*(4.0*res1 + res2)*res1;
477 else if ( i == nx_-1 ) {
478 res1 = (*up)[i-1]-evaluate_target((Real)i*dx_);
479 res2 = (*up)[i]-evaluate_target((Real)(i+1)*dx_);
480 valu += dx_/6.0*(res1 + 4.0*res2)*res2;
483 res1 = (*up)[i-1]-evaluate_target((Real)i*dx_);
484 res2 = (*up)[i]-evaluate_target((Real)(i+1)*dx_);
485 res3 = (*up)[i+1]-evaluate_target((Real)(i+2)*dx_);
486 valu += dx_/6.0*(res1 + 4.0*res2 + res3)*res2;
489 return 0.5*(valu + alpha_*valz);
494 Teuchos::RCP<std::vector<Real> > gup = Teuchos::rcp_const_cast<std::vector<Real> >(
497 Teuchos::RCP<const std::vector<Real> > up =
499 Teuchos::RCP<const std::vector<Real> > zp =
502 std::vector<Real> diff(nx_,0.0);
503 for (
int i=0; i<
nx_; i++) {
504 diff[i] = ((*up)[i]-evaluate_target((Real)(i+1)*dx_));
506 apply_mass(*gup,diff);
511 Teuchos::RCP<std::vector<Real> > gzp = Teuchos::rcp_const_cast<std::vector<Real> >(
514 Teuchos::RCP<const std::vector<Real> > up =
516 Teuchos::RCP<const std::vector<Real> > zp =
519 for (
int i=0; i<nx_+2; i++) {
521 (*gzp)[i] = alpha_*dx_/6.0*(2.0*(*zp)[i]+(*zp)[i+1]);
524 (*gzp)[i] = alpha_*dx_/6.0*(2.0*(*zp)[i]+(*zp)[i-1]);
527 (*gzp)[i] = alpha_*dx_/6.0*((*zp)[i-1]+4.0*(*zp)[i]+(*zp)[i+1]);
534 Teuchos::RCP<std::vector<Real> > hvup = Teuchos::rcp_const_cast<std::vector<Real> >(
537 Teuchos::RCP<const std::vector<Real> > vup =
540 apply_mass(*hvup,*vup);
555 Teuchos::RCP<std::vector<Real> > hvzp = Teuchos::rcp_const_cast<std::vector<Real> >(
558 Teuchos::RCP<const std::vector<Real> > vzp =
561 for (
int i=0; i<nx_+2; i++) {
563 (*hvzp)[i] = alpha_*dx_/6.0*(2.0*(*vzp)[i]+(*vzp)[i+1]);
566 (*hvzp)[i] = alpha_*dx_/6.0*(2.0*(*vzp)[i]+(*vzp)[i-1]);
569 (*hvzp)[i] = alpha_*dx_/6.0*((*vzp)[i-1]+4.0*(*vzp)[i]+(*vzp)[i+1]);
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)
virtual void solve(Vector< Real > &c, Vector< Real > &u, const Vector< Real > &z, Real &tol)
Given , solve for .
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.
Real compute_norm(const std::vector< Real > &r)
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 compute_residual(std::vector< Real > &r, const std::vector< Real > &u, const std::vector< Real > &z)
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 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 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 ...
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 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)
Objective_BurgersControl(Real alpha=1.e-4, int nx=128)
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)
EqualityConstraint_BurgersControl(int nx=128)
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 std::vector< Real > getParameter(void) const
void scale(std::vector< Real > &u, const Real alpha=0.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 apply_mass(std::vector< Real > &Mu, const std::vector< Real > &u)
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 ...