ROL
poisson-control/example_02.cpp
Go to the documentation of this file.
1 
2 // Burgers includes
3 #include "example_02.hpp"
4 // ROL includes
5 #include "ROL_Algorithm.hpp"
6 #include "ROL_StdVector.hpp"
11 #include "ROL_Vector_SimOpt.hpp"
12 // Teuchos includes
13 #include "Teuchos_Time.hpp"
14 #include "Teuchos_oblackholestream.hpp"
15 #include "Teuchos_XMLParameterListHelpers.hpp"
16 #include "Teuchos_GlobalMPISession.hpp"
17 #include "Teuchos_Comm.hpp"
18 #include "Teuchos_DefaultComm.hpp"
19 #include "Teuchos_CommHelpers.hpp"
20 
21 int main( int argc, char *argv[] ) {
22 
23  Teuchos::GlobalMPISession mpiSession(&argc, &argv);
24  Teuchos::RCP<const Teuchos::Comm<int> > comm =
25  Teuchos::DefaultComm<int>::getComm();
26 
27  // This little trick lets us print to std::cout only if a (dummy) command-line argument is provided.
28  int iprint = argc - 1;
29  Teuchos::RCP<std::ostream> outStream;
30  Teuchos::oblackholestream bhs; // outputs nothing
31  if (iprint > 0)
32  outStream = Teuchos::rcp(&std::cout, false);
33  else
34  outStream = Teuchos::rcp(&bhs, false);
35 
36  int errorFlag = 0;
37 
38  // *** Example body.
39 
40  try {
41 
42  /***************************************************************************/
43  /***************** GRAB INPUTS *********************************************/
44  /***************************************************************************/
45  // Get finite element parameter list
46  std::string filename = "example_02.xml";
47  Teuchos::RCP<Teuchos::ParameterList> parlist = Teuchos::rcp( new Teuchos::ParameterList() );
48  Teuchos::updateParametersFromXmlFile( filename, Teuchos::Ptr<Teuchos::ParameterList>(&*parlist) );
49  if ( parlist->get("Display Option",0) && (comm->getRank() > 0) ) {
50  parlist->set("Display Option",0);
51  }
52  // Get ROL parameter list
53  filename = "input.xml";
54  Teuchos::RCP<Teuchos::ParameterList> ROL_parlist = Teuchos::rcp( new Teuchos::ParameterList() );
55  Teuchos::updateParametersFromXmlFile( filename, Teuchos::Ptr<Teuchos::ParameterList>(&*ROL_parlist) );
56 
57  /***************************************************************************/
58  /***************** INITIALIZE SAMPLERS *************************************/
59  /***************************************************************************/
60  int dim = 2;
61  bool useSA = parlist->get("Use Stochastic Approximation",false);
62  int nSamp = 1;
63  if ( !useSA ) {
64  nSamp = parlist->get("Number of Monte Carlo Samples",1000);
65  }
66  std::vector<double> tmp(2); tmp[0] = -1.0; tmp[1] = 1.0;
67  std::vector<std::vector<double> > bounds(dim,tmp);
68  Teuchos::RCP<ROL::BatchManager<double> > bman
69  = Teuchos::rcp(new ROL::StdTeuchosBatchManager<double,int>(comm));
70  Teuchos::RCP<ROL::SampleGenerator<double> > sampler
71  = Teuchos::rcp(new ROL::MonteCarloGenerator<double>(nSamp,bounds,bman,useSA));
72 
73  /***************************************************************************/
74  /***************** INITIALIZE CONTROL VECTOR *******************************/
75  /***************************************************************************/
76  int nx = parlist->get("Number of Elements", 128);
77  Teuchos::RCP<std::vector<double> > z_rcp = Teuchos::rcp( new std::vector<double>(nx+1, 0.0) );
78  Teuchos::RCP<ROL::Vector<double> > z = Teuchos::rcp(new ROL::StdVector<double>(z_rcp));
79  Teuchos::RCP<std::vector<double> > u_rcp = Teuchos::rcp( new std::vector<double>(nx-1, 0.0) );
80  Teuchos::RCP<ROL::Vector<double> > u = Teuchos::rcp(new ROL::StdVector<double>(u_rcp));
82  Teuchos::RCP<std::vector<double> > p_rcp = Teuchos::rcp( new std::vector<double>(nx-1, 0.0) );
83  Teuchos::RCP<ROL::Vector<double> > p = Teuchos::rcp(new ROL::StdVector<double>(p_rcp));
84  Teuchos::RCP<std::vector<double> > U_rcp = Teuchos::rcp( new std::vector<double>(nx+1, 35.0) );
85  Teuchos::RCP<ROL::Vector<double> > U = Teuchos::rcp(new ROL::StdVector<double>(U_rcp));
86  Teuchos::RCP<std::vector<double> > L_rcp = Teuchos::rcp( new std::vector<double>(nx+1, -5.0) );
87  Teuchos::RCP<ROL::Vector<double> > L = Teuchos::rcp(new ROL::StdVector<double>(L_rcp));
89 
90  /***************************************************************************/
91  /***************** INITIALIZE OBJECTIVE FUNCTION ***************************/
92  /***************************************************************************/
93  double alpha = parlist->get("Penalty Parameter", 1.e-4);
94  Teuchos::RCP<FEM<double> > fem = Teuchos::rcp(new FEM<double>(nx));
95  Teuchos::RCP<ROL::Objective_SimOpt<double> > pObj
96  = Teuchos::rcp(new DiffusionObjective<double>(fem, alpha));
97  Teuchos::RCP<ROL::EqualityConstraint_SimOpt<double> > pCon
98  = Teuchos::rcp(new DiffusionEqualityConstraint<double>(fem));
99  Teuchos::RCP<ROL::Objective<double> > robj
100  = Teuchos::rcp(new ROL::Reduced_Objective_SimOpt<double>(pObj,pCon,u,p));
101  ROL::RiskNeutralObjective<double> obj(robj,sampler);
102 
103  /***************************************************************************/
104  /***************** RUN DERIVATIVE CHECK ************************************/
105  /***************************************************************************/
106  if (parlist->get("Run Derivative Check",false)) {
107  // Direction to test finite differences
108  Teuchos::RCP<std::vector<double> > dz_rcp = Teuchos::rcp( new std::vector<double>(nx+1, 0.0) );
109  Teuchos::RCP<ROL::Vector<double> > dz = Teuchos::rcp(new ROL::StdVector<double>(dz_rcp));
110  Teuchos::RCP<std::vector<double> > du_rcp = Teuchos::rcp( new std::vector<double>(nx-1, 0.0) );
111  Teuchos::RCP<ROL::Vector<double> > du = Teuchos::rcp(new ROL::StdVector<double>(du_rcp));
113  // Set to random vectors
114  srand(12345);
115  for (int i=0; i<nx+1; i++) {
116  (*dz_rcp)[i] = 2.0*(double)rand()/(double)RAND_MAX - 1.0;
117  (*z_rcp)[i] = 2.0*(double)rand()/(double)RAND_MAX - 1.0;
118  }
119  for (int i=0; i<nx-1; i++) {
120  (*du_rcp)[i] = 2.0*(double)rand()/(double)RAND_MAX - 1.0;
121  (*u_rcp)[i] = 2.0*(double)rand()/(double)RAND_MAX - 1.0;
122  }
123  // Run derivative checks
124  std::vector<double> param(dim,0.0);
125  robj->setParameter(param);
126  if ( comm->getRank() == 0 ) {
127  std::cout << "\nRUN DERIVATIVE CHECK FOR PARAMETRIZED OBJECTIVE FUNCTION SIMOPT\n";
128  }
129  pObj->checkGradient(x,d,(comm->getRank()==0));
130  pObj->checkHessVec(x,d,(comm->getRank()==0));
131  if ( comm->getRank() == 0 ) {
132  std::cout << "\nRUN DERIVATIVE CHECK FOR PARAMETRIZED EQUALITY CONSTRAINT SIMOPT\n";
133  }
134  pCon->checkApplyJacobian(x,d,*p,(comm->getRank()==0));
135  pCon->checkApplyAdjointJacobian(x,*du,*p,x,(comm->getRank()==0));
136  pCon->checkApplyAdjointHessian(x,*du,d,x,(comm->getRank()==0));
137  if ( comm->getRank() == 0 ) {
138  std::cout << "\nRUN DERIVATIVE CHECK FOR PARAMETRIZED OBJECTIVE FUNCTION\n";
139  }
140  robj->checkGradient(*z,*dz,(comm->getRank()==0));
141  robj->checkHessVec(*z,*dz,(comm->getRank()==0));
142  // Run derivative checks
143  if ( comm->getRank() == 0 ) {
144  std::cout << "\nRUN DERIVATIVE CHECK FOR RISK-NEUTRAL OBJECTIVE FUNCTION\n";
145  }
146  obj.checkGradient(*z,*dz,(comm->getRank()==0));
147  obj.checkHessVec(*z,*dz,(comm->getRank()==0));
148  }
149 
150  /***************************************************************************/
151  /***************** INITIALIZE ROL ALGORITHM ********************************/
152  /***************************************************************************/
153  Teuchos::RCP<ROL::Algorithm<double> > algo;
154  if ( useSA ) {
155  ROL_parlist->sublist("General").set("Recompute Objective Function",false);
156  ROL_parlist->sublist("Step").sublist("Line Search").set("Initial Step Size",0.1/alpha);
157  ROL_parlist->sublist("Step").sublist("Line Search").set("User Defined Initial Step Size",true);
158  ROL_parlist->sublist("Step").sublist("Line Search").sublist("Line-Search Method").set("Type","Iteration Scaling");
159  ROL_parlist->sublist("Step").sublist("Line Search").sublist("Descent Method").set("Type","Steepest Descent");
160  ROL_parlist->sublist("Step").sublist("Line Search").sublist("Curvature Condition").set("Type","Null Curvature Condition");
161  algo = Teuchos::rcp(new ROL::Algorithm<double>("Line Search",*ROL_parlist,false));
162  }
163  else {
164  algo = Teuchos::rcp(new ROL::Algorithm<double>("Trust Region",*ROL_parlist,false));
165  }
166 
167  /***************************************************************************/
168  /***************** PERFORM OPTIMIZATION ************************************/
169  /***************************************************************************/
170  Teuchos::Time timer("Optimization Time",true);
171  z->zero();
172  algo->run(*z,obj,bnd,(comm->getRank()==0));
173  double optTime = timer.stop();
174 
175  /***************************************************************************/
176  /***************** PRINT RESULTS *******************************************/
177  /***************************************************************************/
178  int my_number_samples = sampler->numMySamples(), number_samples = 0;
179  Teuchos::reduceAll<int,int>(*comm,Teuchos::REDUCE_SUM,1,&my_number_samples,&number_samples);
180  int my_number_solves = Teuchos::rcp_dynamic_cast<DiffusionEqualityConstraint<double> >(pCon)->getNumSolves(), number_solves = 0;
181  Teuchos::reduceAll<int,int>(*comm,Teuchos::REDUCE_SUM,1,&my_number_solves,&number_solves);
182  if (comm->getRank() == 0) {
183  std::cout << "Number of Samples = " << number_samples << "\n";
184  std::cout << "Number of Solves = " << number_solves << "\n";
185  std::cout << "Optimization Time = " << optTime << "\n\n";
186  }
187 
188  if ( comm->getRank() == 0 ) {
189  std::ofstream file;
190  if (useSA) {
191  file.open("control_SA.txt");
192  }
193  else {
194  file.open("control_SAA.txt");
195  }
196  std::vector<double> xmesh(fem->nz(),0.0);
197  fem->build_mesh(xmesh);
198  for (int i = 0; i < fem->nz(); i++ ) {
199  file << std::setprecision(std::numeric_limits<double>::digits10) << std::scientific << xmesh[i] << " "
200  << std::setprecision(std::numeric_limits<double>::digits10) << std::scientific << (*z_rcp)[i]
201  << "\n";
202  }
203  file.close();
204  }
205  }
206  catch (std::logic_error err) {
207  *outStream << err.what() << "\n";
208  errorFlag = -1000;
209  }; // end try
210 
211  if (errorFlag != 0)
212  std::cout << "End Result: TEST FAILED\n";
213  else
214  std::cout << "End Result: TEST PASSED\n";
215 
216  return 0;
217 }
218 
219 
220 
221 
Defines the linear algebra or vector space interface for simulation-based optimization.
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.
Provides the std::vector implementation of the ROL::Vector interface.
Provides an interface to run optimization algorithms.
Provides the interface to apply upper and lower bound constraints.
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.
int main(int argc, char *argv[])