ROL
burgers-control/example_02.cpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Rapid Optimization Library (ROL) Package
4 //
5 // Copyright 2014 NTESS and the ROL contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
19 #include "ROL_Stream.hpp"
20 
21 #include "Teuchos_GlobalMPISession.hpp"
22 #include "Teuchos_LAPACK.hpp"
23 
24 #include <iostream>
25 #include <algorithm>
26 
27 #include "example_02.hpp"
28 
29 typedef double RealT;
30 
31 int main(int argc, char *argv[]) {
32 
33  Teuchos::GlobalMPISession mpiSession(&argc, &argv);
34 
35  // This little trick lets us print to std::cout only if a (dummy) command-line argument is provided.
36  int iprint = argc - 1;
37  ROL::Ptr<std::ostream> outStream;
38  ROL::nullstream bhs; // outputs nothing
39  if (iprint > 0)
40  outStream = ROL::makePtrFromRef(std::cout);
41  else
42  outStream = ROL::makePtrFromRef(bhs);
43 
44  int errorFlag = 0;
45 
46  // *** Example body.
47 
48  try {
49  // Initialize full objective function.
50  int nx = 256; // Set spatial discretization.
51  RealT alpha = 1.e-3; // Set penalty parameter.
52  RealT nu = 1e-2; // Viscosity parameter.
53  Objective_BurgersControl<RealT> obj(alpha,nx);
54  // Initialize equality constraints
56  ROL::ParameterList list;
57  list.sublist("SimOpt").sublist("Solve").set("Absolute Residual Tolerance",1.e2*ROL::ROL_EPSILON<RealT>());
58  con.setSolveParameters(list);
59  // Initialize iteration vectors.
60  ROL::Ptr<std::vector<RealT> > z_ptr = ROL::makePtr<std::vector<RealT>>(nx+2, 1.0);
61  ROL::Ptr<std::vector<RealT> > gz_ptr = ROL::makePtr<std::vector<RealT>>(nx+2, 1.0);
62  ROL::Ptr<std::vector<RealT> > yz_ptr = ROL::makePtr<std::vector<RealT>>(nx+2, 1.0);
63  for (int i=0; i<nx+2; i++) {
64  (*z_ptr)[i] = (RealT)rand()/(RealT)RAND_MAX;
65  (*yz_ptr)[i] = (RealT)rand()/(RealT)RAND_MAX;
66  }
67  ROL::StdVector<RealT> z(z_ptr);
68  ROL::StdVector<RealT> gz(gz_ptr);
69  ROL::StdVector<RealT> yz(yz_ptr);
70  ROL::Ptr<ROL::Vector<RealT> > zp = ROL::makePtrFromRef(z);
71  ROL::Ptr<ROL::Vector<RealT> > gzp = ROL::makePtrFromRef(z);
72  ROL::Ptr<ROL::Vector<RealT> > yzp = ROL::makePtrFromRef(yz);
73 
74  ROL::Ptr<std::vector<RealT> > u_ptr = ROL::makePtr<std::vector<RealT>>(nx, 1.0);
75  ROL::Ptr<std::vector<RealT> > gu_ptr = ROL::makePtr<std::vector<RealT>>(nx, 1.0);
76  ROL::Ptr<std::vector<RealT> > yu_ptr = ROL::makePtr<std::vector<RealT>>(nx, 1.0);
77  for (int i=0; i<nx; i++) {
78  (*u_ptr)[i] = (RealT)rand()/(RealT)RAND_MAX;
79  (*yu_ptr)[i] = (RealT)rand()/(RealT)RAND_MAX;
80  }
81  ROL::StdVector<RealT> u(u_ptr);
82  ROL::StdVector<RealT> gu(gu_ptr);
83  ROL::StdVector<RealT> yu(yu_ptr);
84  ROL::Ptr<ROL::Vector<RealT> > up = ROL::makePtrFromRef(u);
85  ROL::Ptr<ROL::Vector<RealT> > gup = ROL::makePtrFromRef(gu);
86  ROL::Ptr<ROL::Vector<RealT> > yup = ROL::makePtrFromRef(yu);
87 
88  ROL::Ptr<std::vector<RealT> > c_ptr = ROL::makePtr<std::vector<RealT>>(nx, 1.0);
89  ROL::Ptr<std::vector<RealT> > l_ptr = ROL::makePtr<std::vector<RealT>>(nx, 1.0);
90  ROL::StdVector<RealT> c(c_ptr);
91  ROL::StdVector<RealT> l(l_ptr);
92 
94  ROL::Vector_SimOpt<RealT> g(gup,gzp);
95  ROL::Vector_SimOpt<RealT> y(yup,yzp);
96 
97  // Check derivatives.
98  obj.checkGradient(x,x,y,true,*outStream);
99  obj.checkHessVec(x,x,y,true,*outStream);
100  con.checkApplyJacobian(x,y,c,true,*outStream);
101  con.checkApplyAdjointJacobian(x,yu,c,x,true,*outStream);
102  con.checkApplyAdjointHessian(x,yu,y,x,true,*outStream);
103 
104  // Initialize reduced objective function.
105  ROL::Ptr<std::vector<RealT> > p_ptr = ROL::makePtr<std::vector<RealT>>(nx, 1.0);
106  ROL::StdVector<RealT> p(p_ptr);
107  ROL::Ptr<ROL::Vector<RealT> > pp = ROL::makePtrFromRef(p);
108  ROL::Ptr<ROL::Objective_SimOpt<RealT> > pobj = ROL::makePtrFromRef(obj);
109  ROL::Ptr<ROL::Constraint_SimOpt<RealT> > pcon = ROL::makePtrFromRef(con);
110  ROL::Reduced_Objective_SimOpt<RealT> robj(pobj,pcon,up,zp,pp);
111  // Check derivatives.
112  robj.checkGradient(z,z,yz,true,*outStream);
113  robj.checkHessVec(z,z,yz,true,*outStream);
114 
115  // Get parameter list.
116  std::string filename = "input.xml";
117  auto parlist = ROL::getParametersFromXmlFile( filename );
118  parlist->sublist("Status Test").set("Gradient Tolerance",1.e-14);
119  parlist->sublist("Status Test").set("Constraint Tolerance",1.e-14);
120  parlist->sublist("Status Test").set("Step Tolerance",1.e-16);
121  parlist->sublist("Status Test").set("Iteration Limit",1000);
122 
123  // Run equality-constrained optimization.
124  RealT zerotol = std::sqrt(ROL::ROL_EPSILON<RealT>());
125  z.zero();
126  con.solve(c,u,z,zerotol);
127  c.zero(); l.zero();
128  {
129  // Define algorithm.
131  // Run Algorithm
132  algo.run(x, obj, con, l, *outStream);
133  }
134  ROL::Ptr<ROL::Vector<RealT> > zCS = z.clone();
135  zCS->set(z);
136 
137  // Run unconstrained optimization.
138  z.zero();
139  {
140  // Define algorithm.
142  // Run Algorithm
143  algo.run(z, z.dual(), robj, *outStream);
144  }
145 
146  // Check solutions.
147  ROL::Ptr<ROL::Vector<RealT> > err = z.clone();
148  err->set(*zCS); err->axpy(-1.,z);
149  errorFlag += ((err->norm()) > 1.e-8) ? 1 : 0;
150  }
151  catch (std::logic_error& err) {
152  *outStream << err.what() << "\n";
153  errorFlag = -1000;
154  }; // end try
155 
156  if (errorFlag != 0)
157  std::cout << "End Result: TEST FAILED\n";
158  else
159  std::cout << "End Result: TEST PASSED\n";
160 
161  return 0;
162 
163 }
164 
virtual const Vector & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
Definition: ROL_Vector.hpp:192
Defines the linear algebra or vector space interface for simulation-based optimization.
virtual std::vector< std::vector< Real > > checkApplyAdjointJacobian(const Vector< Real > &x, const Vector< Real > &v, const Vector< Real > &c, const Vector< Real > &ajv, const bool printToStream=true, std::ostream &outStream=std::cout, const int numSteps=ROL_NUM_CHECKDERIV_STEPS)
Finite-difference check for the application of the adjoint of constraint Jacobian.
void solve(ROL::Vector< Real > &c, ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Given , solve for .
Definition: example_03.hpp:453
virtual void zero()
Set to zero vector.
Definition: ROL_Vector.hpp:133
Defines a no-output stream class ROL::NullStream and a function makeStreamPtr which either wraps a re...
virtual std::vector< std::vector< Real > > checkGradient(const Vector< Real > &x, const Vector< Real > &d, const bool printToStream=true, std::ostream &outStream=std::cout, const int numSteps=ROL_NUM_CHECKDERIV_STEPS, const int order=1)
Finite-difference gradient check.
virtual void run(Vector< Real > &x, const Vector< Real > &g, Objective< Real > &obj, Constraint< Real > &econ, Vector< Real > &emul, const Vector< Real > &eres, std::ostream &outStream=std::cout) override
Run algorithm on equality constrained problems (Type-E). This general interface supports the use of d...
Provides the ROL::Vector interface for scalar values, to be used, for example, with scalar constraint...
Provides an interface to run equality constrained optimization algorithms using the Composite-Step Tr...
void run(Vector< Real > &x, const Vector< Real > &g, Objective< Real > &obj, std::ostream &outStream=std::cout) override
Run algorithm on unconstrained problems (Type-U). This general interface supports the use of dual opt...
virtual Ptr< Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
virtual std::vector< std::vector< Real > > checkApplyAdjointHessian(const Vector< Real > &x, const Vector< Real > &u, const Vector< Real > &v, const Vector< Real > &hv, const std::vector< Real > &step, const bool printToScreen=true, std::ostream &outStream=std::cout, const int order=1)
Finite-difference check for the application of the adjoint of constraint Hessian. ...
virtual std::vector< std::vector< Real > > checkApplyJacobian(const Vector< Real > &x, const Vector< Real > &v, const Vector< Real > &jv, const std::vector< Real > &steps, const bool printToStream=true, std::ostream &outStream=std::cout, const int order=1)
Finite-difference check for the constraint Jacobian application.
basic_nullstream< char, char_traits< char >> nullstream
Definition: ROL_Stream.hpp:38
int main(int argc, char *argv[])
virtual std::vector< std::vector< Real > > checkHessVec(const Vector< Real > &x, const Vector< Real > &v, const bool printToStream=true, std::ostream &outStream=std::cout, const int numSteps=ROL_NUM_CHECKDERIV_STEPS, const int order=1)
Finite-difference Hessian-applied-to-vector check.
virtual void setSolveParameters(ParameterList &parlist)
Set solve parameters.
Provides an interface to run trust-region methods for unconstrained optimization algorithms.