ROL
example_04.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 
15 #include "ROL_Algorithm.hpp"
18 #include "ROL_Vector_SimOpt.hpp"
19 #include "ROL_ParameterList.hpp"
20 
21 #include "ROL_Stream.hpp"
22 #include "Teuchos_GlobalMPISession.hpp"
23 
24 #include <iostream>
25 #include <algorithm>
26 
27 #include "example_04.hpp"
28 
29 typedef double RealT;
36 
37 int main(int argc, char *argv[]) {
38 
39  Teuchos::GlobalMPISession mpiSession(&argc, &argv);
40  // This little trick lets us print to std::cout only if a (dummy) command-line argument is provided.
41  int iprint = argc - 1;
42  ROL::Ptr<std::ostream> outStream;
43  ROL::nullstream bhs; // outputs nothing
44  if (iprint > 0)
45  outStream = ROL::makePtrFromRef(std::cout);
46  else
47  outStream = ROL::makePtrFromRef(bhs);
48 
49  int errorFlag = 0;
50 
51  // *** Example body.
52  try {
53  /*************************************************************************/
54  /************* INITIALIZE BURGERS FEM CLASS ******************************/
55  /*************************************************************************/
56  int nx = 128; // Set spatial discretization.
57  RealT alpha = 1.e-3; // Set penalty parameter.
58  RealT nu = 1e-2; // Viscosity parameter.
59  RealT nl = 1.0; // Nonlinearity parameter (1 = Burgers, 0 = linear).
60  RealT u0 = 1.0; // Dirichlet boundary condition at x=0.
61  RealT u1 = 0.0; // Dirichlet boundary condition at x=1.
62  RealT f = 0.0; // Constant volumetric force.
63  RealT cH1 = 1.0; // Scale for derivative term in H1 norm.
64  RealT cL2 = 0.0; // Scale for mass term in H1 norm.
65  ROL::Ptr<BurgersFEM<RealT> > fem
66  = ROL::makePtr<BurgersFEM<RealT>>(nx,nu,nl,u0,u1,f,cH1,cL2);
67  fem->test_inverse_mass(*outStream);
68  fem->test_inverse_H1(*outStream);
69  /*************************************************************************/
70  /************* INITIALIZE SIMOPT OBJECTIVE FUNCTION **********************/
71  /*************************************************************************/
72  ROL::Ptr<std::vector<RealT> > ud_ptr
73  = ROL::makePtr<std::vector<RealT>>(nx, 1.);
74  ROL::Ptr<ROL::Vector<RealT> > ud
75  = ROL::makePtr<L2VectorPrimal<RealT>>(ud_ptr,fem);
76  Objective_BurgersControl<RealT> obj(fem,ud,alpha);
77  /*************************************************************************/
78  /************* INITIALIZE SIMOPT EQUALITY CONSTRAINT *********************/
79  /*************************************************************************/
80  bool useEChessian = true;
81  Constraint_BurgersControl<RealT> con(fem, useEChessian);
82  /*************************************************************************/
83  /************* INITIALIZE BOUND CONSTRAINTS ******************************/
84  /*************************************************************************/
85  // INITIALIZE STATE CONSTRAINTS
86  std::vector<RealT> Ulo(nx, 0.), Uhi(nx, 1.);
87  //std::vector<RealT> Ulo(nx, -1.e8), Uhi(nx, 1.e8);
88  ROL::Ptr<ROL::BoundConstraint<RealT> > Ubnd
89  = ROL::makePtr<H1BoundConstraint<RealT>>(Ulo,Uhi,fem);
90  //Ubnd->deactivate();
91  // INITIALIZE CONTROL CONSTRAINTS
92  //std::vector<RealT> Zlo(nx+2, -1.e8), Zhi(nx+2, 1.e8);
93  std::vector<RealT> Zlo(nx+2,0.), Zhi(nx+2,2.);
94  ROL::Ptr<ROL::BoundConstraint<RealT> > Zbnd
95  = ROL::makePtr<L2BoundConstraint<RealT>>(Zlo,Zhi,fem);
96  //Zbnd->deactivate();
97  // INITIALIZE SIMOPT BOUND CONSTRAINTS
99  bnd.deactivate();
100  /*************************************************************************/
101  /************* INITIALIZE VECTOR STORAGE *********************************/
102  /*************************************************************************/
103  // INITIALIZE CONTROL VECTORS
104  ROL::Ptr<std::vector<RealT> > z_ptr
105  = ROL::makePtr<std::vector<RealT>>(nx+2, 0.);
106  ROL::Ptr<std::vector<RealT> > zrand_ptr
107  = ROL::makePtr<std::vector<RealT>>(nx+2, 1.);
108  ROL::Ptr<std::vector<RealT> > gz_ptr
109  = ROL::makePtr<std::vector<RealT>>(nx+2, 1.);
110  ROL::Ptr<std::vector<RealT> > yz_ptr
111  = ROL::makePtr<std::vector<RealT>>(nx+2, 1.);
112  for (int i=0; i<nx+2; i++) {
113  (*zrand_ptr)[i] = 10.*(RealT)rand()/(RealT)RAND_MAX-5.;
114  (*yz_ptr)[i] = 10.*(RealT)rand()/(RealT)RAND_MAX-5.;
115  }
116  ROL::Ptr<ROL::Vector<RealT> > zp
117  = ROL::makePtr<PrimalControlVector>(z_ptr,fem);
118  ROL::Ptr<ROL::Vector<RealT> > zrandp
119  = ROL::makePtr<PrimalControlVector>(zrand_ptr,fem);
120  ROL::Ptr<ROL::Vector<RealT> > gzp
121  = ROL::makePtr<DualControlVector>(gz_ptr,fem);
122  ROL::Ptr<ROL::Vector<RealT> > yzp
123  = ROL::makePtr<PrimalControlVector>(yz_ptr,fem);
124  // INITIALIZE STATE VECTORS
125  ROL::Ptr<std::vector<RealT> > u_ptr
126  = ROL::makePtr<std::vector<RealT>>(nx, 1.);
127  ROL::Ptr<std::vector<RealT> > gu_ptr
128  = ROL::makePtr<std::vector<RealT>>(nx, 1.);
129  ROL::Ptr<std::vector<RealT> > yu_ptr
130  = ROL::makePtr<std::vector<RealT>>(nx, 1.);
131  for (int i=0; i<nx; i++) {
132  (*yu_ptr)[i] = 10.*(RealT)rand()/(RealT)RAND_MAX-5.;
133  }
134  ROL::Ptr<ROL::Vector<RealT> > up
135  = ROL::makePtr<PrimalStateVector>(u_ptr,fem);
136  ROL::Ptr<ROL::Vector<RealT> > gup
137  = ROL::makePtr<DualStateVector>(gu_ptr,fem);
138  ROL::Ptr<ROL::Vector<RealT> > yup
139  = ROL::makePtr<PrimalStateVector>(yu_ptr,fem);
140  // INITIALIZE CONSTRAINT VECTORS
141  ROL::Ptr<std::vector<RealT> > c_ptr
142  = ROL::makePtr<std::vector<RealT>>(nx, 1.);
143  ROL::Ptr<std::vector<RealT> > l_ptr
144  = ROL::makePtr<std::vector<RealT>>(nx, 1.);
145  for (int i=0; i<nx; i++) {
146  (*l_ptr)[i] = (RealT)rand()/(RealT)RAND_MAX;
147  }
148  PrimalConstraintVector c(c_ptr,fem);
149  DualConstraintVector l(l_ptr,fem);
150  // INITIALIZE SIMOPT VECTORS
151  ROL::Vector_SimOpt<RealT> x(up,zp);
152  ROL::Vector_SimOpt<RealT> g(gup,gzp);
153  ROL::Vector_SimOpt<RealT> y(yup,yzp);
154  // READ IN XML INPUT
155  std::string filename = "input.xml";
156  auto parlist = ROL::getParametersFromXmlFile( filename );
157 
158  /*************************************************************************/
159  /************* CHECK DERIVATIVES AND CONSISTENCY *************************/
160  /*************************************************************************/
161  zp->set(*zrandp);
162  // CHECK OBJECTIVE DERIVATIVES
163  obj.checkGradient(x,g,y,true,*outStream);
164  obj.checkHessVec(x,g,y,true,*outStream);
165  // CHECK EQUALITY CONSTRAINT DERIVATIVES
166  con.checkApplyJacobian(x,y,c,true,*outStream);
167  con.checkApplyAdjointHessian(x,*yup,y,g,true,*outStream);
168  // CHECK EQUALITY CONSTRAINT CONSISTENCY
169  con.checkSolve(*up,*zp,c,true,*outStream);
170  con.checkAdjointConsistencyJacobian_1(l,*yup,*up,*zp,true,*outStream);
171  con.checkAdjointConsistencyJacobian_2(l,*yzp,*up,*zp,true,*outStream);
172  con.checkInverseJacobian_1(c,*yup,*up,*zp,true,*outStream);
173  con.checkInverseAdjointJacobian_1(c,*yup,*up,*zp,true,*outStream);
174  *outStream << "\n";
175  // CHECK PENALTY OBJECTIVE DERIVATIVES
176  ROL::Ptr<ROL::Objective<RealT> > obj_ptr = ROL::makePtrFromRef(obj);
177  ROL::Ptr<ROL::Constraint<RealT> > con_ptr = ROL::makePtrFromRef(con);
178  ROL::Ptr<ROL::BoundConstraint<RealT> > bnd_ptr = ROL::makePtrFromRef(bnd);
179  ROL::MoreauYosidaPenalty<RealT> myPen(obj_ptr,bnd_ptr,x,*parlist);
180  myPen.checkGradient(x, y, true, *outStream);
181  myPen.checkHessVec(x, g, y, true, *outStream);
182  ROL::AugmentedLagrangian<RealT> myAugLag(obj_ptr,con_ptr,l,1.,x,c,*parlist);
183  myAugLag.checkGradient(x, y, true, *outStream);
184  myAugLag.checkHessVec(x, g, y, true, *outStream);
185  /*************************************************************************/
186  /************* RUN OPTIMIZATION ******************************************/
187  /*************************************************************************/
188  // SOLVE USING MOREAU-YOSIDA PENALTY
189  ROL::Ptr<ROL::Step<RealT>>
190  stepMY = ROL::makePtr<ROL::MoreauYosidaPenaltyStep<RealT>>(*parlist);
191  ROL::Ptr<ROL::StatusTest<RealT>>
192  statusMY = ROL::makePtr<ROL::ConstraintStatusTest<RealT>>(*parlist);
193  ROL::Algorithm<RealT> algoMY(stepMY,statusMY,false);
194  zp->set(*zrandp);
195  RealT zerotol = std::sqrt(ROL::ROL_EPSILON<RealT>());
196  con.solve(c,*up,*zp,zerotol);
197  obj.gradient_1(*gup,*up,*zp,zerotol);
198  gup->scale(-1.0);
199  con.applyInverseAdjointJacobian_1(l,*gup,*up,*zp,zerotol);
200  gup->zero(); c.zero();
201  algoMY.run(x, g, l, c, myPen, con, bnd, true, *outStream);
202  ROL::Ptr<ROL::Vector<RealT> > xMY = x.clone();
203  xMY->set(x);
204  // SOLVE USING AUGMENTED LAGRANGIAN
205  ROL::Ptr<ROL::Step<RealT>>
206  stepAL = ROL::makePtr<ROL::AugmentedLagrangianStep<RealT>>(*parlist);
207  ROL::Ptr<ROL::StatusTest<RealT>>
208  statusAL = ROL::makePtr<ROL::ConstraintStatusTest<RealT>>(*parlist);
209  ROL::Algorithm<RealT> algoAL(stepAL,statusAL,false);
210  zp->set(*zrandp);
211  con.solve(c,*up,*zp,zerotol);
212  obj.gradient_1(*gup,*up,*zp,zerotol);
213  gup->scale(-1.0);
214  con.applyInverseAdjointJacobian_1(l,*gup,*up,*zp,zerotol);
215  gup->zero(); c.zero();
216  algoAL.run(x, g, l, c, myAugLag, con, bnd, true, *outStream);
217  // COMPARE SOLUTIONS
218  ROL::Ptr<ROL::Vector<RealT> > err = x.clone();
219  err->set(x); err->axpy(-1.,*xMY);
220  errorFlag += ((err->norm() > 1.e-7*x.norm()) ? 1 : 0);
221  }
222  catch (std::logic_error& err) {
223  *outStream << err.what() << "\n";
224  errorFlag = -1000;
225  }; // end try
226 
227  if (errorFlag != 0)
228  std::cout << "End Result: TEST FAILED\n";
229  else
230  std::cout << "End Result: TEST PASSED\n";
231 
232  return 0;
233 }
Provides the interface to evaluate the augmented Lagrangian.
virtual Real checkSolve(const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &c, const bool printToStream=true, std::ostream &outStream=std::cout)
Defines the linear algebra or vector space interface for simulation-based optimization.
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 std::vector< std::string > run(Vector< Real > &x, Objective< Real > &obj, bool print=false, std::ostream &outStream=std::cout, bool printVectors=false, std::ostream &vectorStream=std::cout)
Run algorithm on unconstrained problems (Type-U). This is the primary Type-U interface.
virtual Real checkAdjointConsistencyJacobian_2(const Vector< Real > &w, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, const bool printToStream=true, std::ostream &outStream=std::cout)
Check the consistency of the Jacobian and its adjoint. This is the primary interface.
L2VectorPrimal< RealT > PrimalControlVector
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 ...
Definition: test_04.hpp:772
Real norm() const
Returns where .
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.
ROL::Ptr< Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
virtual Real checkInverseJacobian_1(const Vector< Real > &jv, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, const bool printToStream=true, std::ostream &outStream=std::cout)
H1VectorDual< RealT > DualStateVector
L2VectorDual< RealT > DualControlVector
Provides an interface to run optimization algorithms.
Provides the interface to evaluate the Moreau-Yosida penalty function.
H1VectorDual< RealT > PrimalConstraintVector
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[])
Provides definitions of equality constraint and objective for example_04.
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 Real checkAdjointConsistencyJacobian_1(const Vector< Real > &w, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, const bool printToStream=true, std::ostream &outStream=std::cout)
Check the consistency of the Jacobian and its adjoint. This is the primary interface.
H1VectorPrimal< RealT > PrimalStateVector
virtual Real checkInverseAdjointJacobian_1(const Vector< Real > &jv, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, const bool printToStream=true, std::ostream &outStream=std::cout)
H1VectorPrimal< RealT > DualConstraintVector