11 #include "Teuchos_LAPACK.hpp"
12 #include "Teuchos_GlobalMPISession.hpp"
13 #include "Teuchos_Comm.hpp"
14 #include "Teuchos_DefaultComm.hpp"
15 #include "Teuchos_CommHelpers.hpp"
17 #include "ROL_ParameterList.hpp"
29 #include "ROL_StdTeuchosBatchManager.hpp"
38 return std::sqrt(
dot(r,r));
41 Real
dot(
const std::vector<Real> &x,
const std::vector<Real> &y) {
43 Real c = (((int)x.size()==
nx_) ? 4.0 : 2.0);
44 for (
unsigned i=0; i<x.size(); i++) {
46 ip +=
dx_/6.0*(c*x[i] + x[i+1])*y[i];
48 else if ( i == x.size()-1 ) {
49 ip +=
dx_/6.0*(x[i-1] + c*x[i])*y[i];
52 ip +=
dx_/6.0*(x[i-1] + 4.0*x[i] + x[i+1])*y[i];
60 void update(std::vector<Real> &u,
const std::vector<Real> &s,
const Real alpha=1.0) {
61 for (
unsigned i=0; i<u.size(); i++) {
66 void scale(std::vector<Real> &u,
const Real alpha=0.0) {
67 for (
unsigned i=0; i<u.size(); i++) {
73 const std::vector<Real> &z) {
74 r.clear(); r.resize(
nx_,0.0);
75 const std::vector<Real> param =
77 Real nu = std::pow(10.0,param[0]-2.0);
78 Real f = param[1]/100.0;
79 Real u0 = 1.0+param[2]/1000.0;
80 Real u1 = param[3]/1000.0;
81 for (
int i=0; i<
nx_; i++) {
84 r[i] = nu/
dx_*(2.0*u[i]-u[i+1]);
87 r[i] = nu/
dx_*(2.0*u[i]-u[i-1]);
90 r[i] = nu/
dx_*(2.0*u[i]-u[i-1]-u[i+1]);
94 r[i] += u[i+1]*(u[i]+u[i+1])/6.0;
97 r[i] -= u[i-1]*(u[i-1]+u[i])/6.0;
100 r[i] -=
dx_/6.0*(z[i]+4.0*z[i+1]+z[i+2]);
105 r[ 0] -= u0*u[ 0]/6.0 + u0*u0/6.0 + nu*u0/
dx_;
106 r[nx_-1] += u1*u[nx_-1]/6.0 + u1*u1/6.0 - nu*u1/
dx_;
110 const std::vector<Real> &u) {
111 const std::vector<Real> param =
113 Real nu = std::pow(10.0,param[0]-2.0);
114 Real u0 = 1.0+param[2]/1000.0;
115 Real u1 = param[3]/1000.0;
117 d.clear(); d.resize(
nx_,nu*2.0/
dx_);
118 dl.clear(); dl.resize(
nx_-1,-nu/
dx_);
119 du.clear(); du.resize(
nx_-1,-nu/
dx_);
121 for (
int i=0; i<
nx_; i++) {
123 dl[i] += (-2.0*u[i]-u[i+1])/6.0;
128 du[i-1] += (u[i-1]+2.0*u[i])/6.0;
136 void linear_solve(std::vector<Real> &u, std::vector<Real> &dl, std::vector<Real> &d, std::vector<Real> &du,
137 const std::vector<Real> &r,
const bool transpose =
false) {
138 u.assign(r.begin(),r.end());
140 Teuchos::LAPACK<int,Real> lp;
141 std::vector<Real> du2(
nx_-2,0.0);
142 std::vector<int> ipiv(
nx_,0);
146 lp.GTTRF(
nx_,&dl[0],&d[0],&du[0],&du2[0],&ipiv[0],&info);
151 lp.GTTRS(trans,
nx_,nhrs,&dl[0],&d[0],&du[0],&du2[0],&ipiv[0],&u[0],ldb,&info);
160 ROL::Ptr<std::vector<Real> > cp =
162 ROL::Ptr<const std::vector<Real> > up =
164 ROL::Ptr<const std::vector<Real> > zp =
170 ROL::Ptr<std::vector<Real> > up =
172 up->assign(up->size(),
static_cast<Real
>(1));
178 ROL::Ptr<std::vector<Real> > jvp =
180 ROL::Ptr<const std::vector<Real> > vp =
182 ROL::Ptr<const std::vector<Real> > up =
184 ROL::Ptr<const std::vector<Real> > zp =
186 const std::vector<Real> param =
188 Real nu = std::pow(10.0,param[0]-2.0);
189 Real u0 = 1.0+param[2]/1000.0;
190 Real u1 = param[3]/1000.0;
192 for (
int i = 0; i <
nx_; i++) {
193 (*jvp)[i] = nu/
dx_*2.0*(*vp)[i];
195 (*jvp)[i] += -nu/
dx_*(*vp)[i-1]
196 -(*up)[i-1]/6.0*(*vp)[i]
197 -((*up)[i]+2.0*(*up)[i-1])/6.0*(*vp)[i-1];
200 (*jvp)[i] += -nu/
dx_*(*vp)[i+1]
201 +(*up)[i+1]/6.0*(*vp)[i]
202 +((*up)[i]+2.0*(*up)[i+1])/6.0*(*vp)[i+1];
205 (*jvp)[ 0] -= u0/6.0*(*vp)[0];
206 (*jvp)[nx_-1] += u1/6.0*(*vp)[nx_-1];
211 ROL::Ptr<std::vector<Real> > jvp =
213 ROL::Ptr<const std::vector<Real>> vp =
215 ROL::Ptr<const std::vector<Real>> up =
217 ROL::Ptr<const std::vector<Real>> zp =
219 for (
int i=0; i<
nx_; i++) {
221 (*jvp)[i] = -
dx_/6.0*((*vp)[i]+4.0*(*vp)[i+1]+(*vp)[i+2]);
227 ROL::Ptr<std::vector<Real> > ijvp =
229 ROL::Ptr<const std::vector<Real> > vp =
231 ROL::Ptr<const std::vector<Real> > up =
233 ROL::Ptr<const std::vector<Real> > zp =
236 std::vector<Real> d(
nx_,0.0);
237 std::vector<Real> dl(
nx_-1,0.0);
238 std::vector<Real> du(
nx_-1,0.0);
246 ROL::Ptr<std::vector<Real> > jvp =
248 ROL::Ptr<const std::vector<Real> > vp =
250 ROL::Ptr<const std::vector<Real> > up =
252 ROL::Ptr<const std::vector<Real> > zp =
254 const std::vector<Real> param =
256 Real nu = std::pow(10.0,param[0]-2.0);
257 Real u0 = 1.0+param[2]/1000.0;
258 Real u1 = param[3]/1000.0;
260 for (
int i = 0; i <
nx_; i++) {
261 (*jvp)[i] = nu/
dx_*2.0*(*vp)[i];
263 (*jvp)[i] += -nu/
dx_*(*vp)[i-1]
264 -(*up)[i-1]/6.0*(*vp)[i]
265 +((*up)[i-1]+2.0*(*up)[i])/6.0*(*vp)[i-1];
268 (*jvp)[i] += -nu/
dx_*(*vp)[i+1]
269 +(*up)[i+1]/6.0*(*vp)[i]
270 -((*up)[i+1]+2.0*(*up)[i])/6.0*(*vp)[i+1];
273 (*jvp)[ 0] -= u0/6.0*(*vp)[0];
274 (*jvp)[nx_-1] += u1/6.0*(*vp)[nx_-1];
279 ROL::Ptr<std::vector<Real> > jvp =
281 ROL::Ptr<const std::vector<Real> > vp =
283 ROL::Ptr<const std::vector<Real> > up =
285 ROL::Ptr<const std::vector<Real> > zp =
287 for (
int i=0; i<
nx_+2; i++) {
289 (*jvp)[i] = -
dx_/6.0*(*vp)[i];
292 (*jvp)[i] = -
dx_/6.0*(4.0*(*vp)[i-1]+(*vp)[i]);
294 else if ( i == nx_ ) {
295 (*jvp)[i] = -
dx_/6.0*(4.0*(*vp)[i-1]+(*vp)[i-2]);
297 else if ( i == nx_+1 ) {
298 (*jvp)[i] = -
dx_/6.0*(*vp)[i-2];
301 (*jvp)[i] = -
dx_/6.0*((*vp)[i-2]+4.0*(*vp)[i-1]+(*vp)[i]);
308 ROL::Ptr<std::vector<Real> > iajvp =
310 ROL::Ptr<const std::vector<Real> > vp =
312 ROL::Ptr<const std::vector<Real>> up =
315 std::vector<Real> d(
nx_,0.0);
316 std::vector<Real> du(
nx_-1,0.0);
317 std::vector<Real> dl(
nx_-1,0.0);
325 ROL::Ptr<std::vector<Real> > ahwvp =
327 ROL::Ptr<const std::vector<Real> > wp =
329 ROL::Ptr<const std::vector<Real> > vp =
331 ROL::Ptr<const std::vector<Real> > up =
333 ROL::Ptr<const std::vector<Real> > zp =
335 for (
int i=0; i<
nx_; i++) {
339 (*ahwvp)[i] += ((*wp)[i]*(*vp)[i+1] - (*wp)[i+1]*(2.0*(*vp)[i]+(*vp)[i+1]))/6.0;
342 (*ahwvp)[i] += ((*wp)[i-1]*((*vp)[i-1]+2.0*(*vp)[i]) - (*wp)[i]*(*vp)[i-1])/6.0;
376 case 1: val = ((x<0.5) ? 1.0 : 0.0);
break;
377 case 2: val = 1.0;
break;
378 case 3: val = std::abs(std::sin(8.0*M_PI*x));
break;
379 case 4: val = std::exp(-0.5*(x-0.5)*(x-0.5));
break;
384 Real
dot(
const std::vector<Real> &x,
const std::vector<Real> &y) {
386 Real c = (((int)x.size()==
nx_) ? 4.0 : 2.0);
387 for (
unsigned i=0; i<x.size(); i++) {
389 ip +=
dx_/6.0*(c*x[i] + x[i+1])*y[i];
391 else if ( i == x.size()-1 ) {
392 ip +=
dx_/6.0*(x[i-1] + c*x[i])*y[i];
395 ip +=
dx_/6.0*(x[i-1] + 4.0*x[i] + x[i+1])*y[i];
401 void apply_mass(std::vector<Real> &Mu,
const std::vector<Real> &u ) {
402 Mu.resize(u.size(),0.0);
403 Real c = (((int)u.size()==
nx_) ? 4.0 : 2.0);
404 for (
unsigned i=0; i<u.size(); i++) {
406 Mu[i] =
dx_/6.0*(c*u[i] + u[i+1]);
408 else if ( i == u.size()-1 ) {
409 Mu[i] =
dx_/6.0*(u[i-1] + c*u[i]);
412 Mu[i] =
dx_/6.0*(u[i-1] + 4.0*u[i] + u[i+1]);
423 dx_ = 1.0/((Real)nx+1.0);
429 ROL::Ptr<const std::vector<Real> > up =
431 ROL::Ptr<const std::vector<Real> > zp =
434 Real res1 = 0.0, res2 = 0.0, res3 = 0.0;
435 Real valu = 0.0, valz =
dot(*zp,*zp);
436 for (
int i=0; i<
nx_; i++) {
440 valu +=
dx_/6.0*(4.0*res1 + res2)*res1;
442 else if ( i == nx_-1 ) {
445 valu += dx_/6.0*(res1 + 4.0*res2)*res2;
451 valu += dx_/6.0*(res1 + 4.0*res2 + res3)*res2;
454 return 0.5*(valu +
alpha_*valz);
459 ROL::Ptr<std::vector<Real> > gup =
462 ROL::Ptr<const std::vector<Real> > up =
464 ROL::Ptr<const std::vector<Real> > zp =
467 std::vector<Real> diff(
nx_,0.0);
468 for (
int i=0; i<
nx_; i++) {
476 ROL::Ptr<std::vector<Real> > gzp =
479 ROL::Ptr<const std::vector<Real> > up =
481 ROL::Ptr<const std::vector<Real> > zp =
484 for (
int i=0; i<
nx_+2; i++) {
486 (*gzp)[i] =
alpha_*
dx_/6.0*(2.0*(*zp)[i]+(*zp)[i+1]);
489 (*gzp)[i] =
alpha_*
dx_/6.0*(2.0*(*zp)[i]+(*zp)[i-1]);
492 (*gzp)[i] =
alpha_*
dx_/6.0*((*zp)[i-1]+4.0*(*zp)[i]+(*zp)[i+1]);
499 ROL::Ptr<std::vector<Real> > hvup =
502 ROL::Ptr<const std::vector<Real> > vup =
520 ROL::Ptr<std::vector<Real> > hvzp =
523 ROL::Ptr<const std::vector<Real> > vzp =
526 for (
int i=0; i<
nx_+2; i++) {
528 (*hvzp)[i] =
alpha_*
dx_/6.0*(2.0*(*vzp)[i]+(*vzp)[i+1]);
531 (*hvzp)[i] =
alpha_*
dx_/6.0*(2.0*(*vzp)[i]+(*vzp)[i-1]);
534 (*hvzp)[i] =
alpha_*
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 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 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 ...
const std::vector< Real > getParameter(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 ...
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 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)
Constraint_BurgersControl(int nx=128)
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)
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 .
virtual void solve(Vector< Real > &c, Vector< Real > &u, const Vector< Real > &z, Real &tol)
Given , solve for .
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)