ROL
poisson-control/example_03.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 
10 // Burgers includes
11 #include "example_02.hpp"
12 // ROL includes
13 #include "ROL_StdVector.hpp"
14 #include "ROL_StdTeuchosBatchManager.hpp"
17 #include "ROL_ParameterList.hpp"
18 #include "ROL_Solver.hpp"
20 #include "ROL_PrimalDualRisk.hpp"
21 
22 // Teuchos includes
23 #include "Teuchos_Time.hpp"
24 #include "ROL_Stream.hpp"
25 #include "Teuchos_GlobalMPISession.hpp"
26 #include "Teuchos_Comm.hpp"
27 #include "Teuchos_DefaultComm.hpp"
28 #include "Teuchos_CommHelpers.hpp"
29 
30 int main( int argc, char *argv[] ) {
31 
32  Teuchos::GlobalMPISession mpiSession(&argc, &argv);
33 
34  auto comm = ROL::toPtr( Teuchos::DefaultComm<int>::getComm() );
35 
36  // This little trick lets us print to std::cout only if a (dummy) command-line argument is provided.
37  int iprint = argc - 1;
38  ROL::Ptr<std::ostream> outStream;
39  ROL::nullstream bhs; // outputs nothing
40  if (iprint > 0 && comm->getRank()==0)
41  outStream = ROL::makePtrFromRef(std::cout);
42  else
43  outStream = ROL::makePtrFromRef(bhs);
44 
45  int errorFlag = 0;
46 
47  // *** Example body.
48 
49  try {
50 
51  /***************************************************************************/
52  /***************** GRAB INPUTS *********************************************/
53  /***************************************************************************/
54  // Get finite element parameter list
55  std::string filename = "input_ex03.xml";
56  auto parlist = ROL::getParametersFromXmlFile( filename );
57 
58  if ( parlist->sublist("Problem Data").get("Display Option",0) && (comm->getRank() > 0) ) {
59  parlist->set("Display Option",0);
60  }
61 
62  /***************************************************************************/
63  /***************** INITIALIZE SAMPLERS *************************************/
64  /***************************************************************************/
65  int dim = 2;
66  int nSamp = parlist->sublist("Problem Data").get("Number of Monte Carlo Samples",1000);
67  std::vector<double> tmp(2); tmp[0] = -1.0; tmp[1] = 1.0;
68  std::vector<std::vector<double> > bounds(dim,tmp);
69  ROL::Ptr<ROL::BatchManager<double> > bman
70  = ROL::makePtr<ROL::StdTeuchosBatchManager<double,int>>(comm);
71  ROL::Ptr<ROL::SampleGenerator<double> > sampler
72  = ROL::makePtr<ROL::MonteCarloGenerator<double>>(nSamp,bounds,bman,false);
73 
74  /***************************************************************************/
75  /***************** INITIALIZE CONTROL VECTOR *******************************/
76  /***************************************************************************/
77  int nx = parlist->sublist("Problem Data").get("Number of Elements", 128);
78  ROL::Ptr<std::vector<double> > z_ptr = ROL::makePtr<std::vector<double>>(nx+1, 0.0);
79  ROL::Ptr<ROL::Vector<double> > z = ROL::makePtr<ROL::StdVector<double>>(z_ptr);
80  ROL::Ptr<ROL::Vector<double> > u = ROL::makePtr<ROL::StdVector<double>>(nx-1);
81  ROL::Ptr<ROL::Vector<double> > p = ROL::makePtr<ROL::StdVector<double>>(nx-1);
82 
83  /***************************************************************************/
84  /***************** INITIALIZE OBJECTIVE FUNCTION ***************************/
85  /***************************************************************************/
86  double alpha = parlist->sublist("Problem Data").get("Penalty Parameter", 1.e-4);
87  ROL::Ptr<FEM<double> > fem = ROL::makePtr<FEM<double>>(nx);
88  ROL::Ptr<ROL::Objective_SimOpt<double> > pObj
89  = ROL::makePtr<DiffusionObjective<double>>(fem, alpha);
90  ROL::Ptr<ROL::Constraint_SimOpt<double> > pCon
91  = ROL::makePtr<DiffusionConstraint<double>>(fem);
92  ROL::Ptr<ROL::Objective<double> > robj
93  = ROL::makePtr<ROL::Reduced_Objective_SimOpt<double>>(pObj,pCon,u,z,p);
94  robj->setParameter({0.0,0.0});
95 
96  /***************************************************************************/
97  /***************** INITIALIZE ROL ALGORITHM ********************************/
98  /***************************************************************************/
99  bool runBundle = parlist->sublist("Problem Data").get("Run Bundle",false);
100  // Solve using bundle
101  if (runBundle) {
102  z->zero();
103  ROL::Ptr<ROL::StochasticProblem<double>> problem2
104  = ROL::makePtr<ROL::StochasticProblem<double>>(robj, z);
105  problem2->makeObjectiveStochastic(*parlist, sampler);
106  problem2->finalize(false,true,*outStream);
107  parlist->sublist("Step").set("Type","Bundle");
108  parlist->sublist("Step").sublist("Bundle").set("Distance Measure Coefficient",0.0);
109  ROL::Solver<double> solver2(problem2,*parlist);
110  solver2.solve(*outStream);
111  }
112 
113  ROL::Ptr<ROL::Problem<double>> problem
114  = ROL::makePtr<ROL::Problem<double>>(robj, z);
115  ROL::PrimalDualRisk<double> solver(problem, sampler, *parlist);
116  if (parlist->sublist("Problem Data").get("Run Derivative Check",false)) {
117  problem->check(true,*outStream);
118  solver.check(*outStream);
119  }
120  solver.run(*outStream);
121 
122  /***************************************************************************/
123  /***************** PRINT RESULTS *******************************************/
124  /***************************************************************************/
125  int my_number_samples = sampler->numMySamples(), number_samples = 0;
126  Teuchos::reduceAll<int,int>(*comm,Teuchos::REDUCE_SUM,1,&my_number_samples,&number_samples);
127  int my_number_solves = ROL::dynamicPtrCast<DiffusionConstraint<double> >(pCon)->getNumSolves(), number_solves = 0;
128  Teuchos::reduceAll<int,int>(*comm,Teuchos::REDUCE_SUM,1,&my_number_solves,&number_solves);
129  if (comm->getRank() == 0) {
130  std::cout << "Number of Samples = " << number_samples << "\n";
131  std::cout << "Number of Solves = " << number_solves << "\n";
132  }
133 
134  if ( comm->getRank() == 0 ) {
135  std::ofstream file;
136  file.open("control.txt");
137  std::vector<double> xmesh(fem->nz(),0.0);
138  fem->build_mesh(xmesh);
139  for (int i = 0; i < fem->nz(); i++ ) {
140  file << std::setprecision(std::numeric_limits<double>::digits10) << std::scientific << xmesh[i] << " "
141  << std::setprecision(std::numeric_limits<double>::digits10) << std::scientific << (*z_ptr)[i]
142  << "\n";
143  }
144  file.close();
145  }
146  }
147  catch (std::logic_error& err) {
148  *outStream << err.what() << "\n";
149  errorFlag = -1000;
150  }; // end try
151 
152  if (errorFlag != 0)
153  std::cout << "End Result: TEST FAILED\n";
154  else
155  std::cout << "End Result: TEST PASSED\n";
156 
157  return 0;
158 }
159 
160 
161 
162 
Provides a simplified interface for solving a wide range of optimization problems.
Definition: ROL_Solver.hpp:30
Defines a no-output stream class ROL::NullStream and a function makeStreamPtr which either wraps a re...
int solve(const Ptr< StatusTest< Real >> &status=nullPtr, bool combineStatus=true)
Solve optimization problem with no iteration output.
basic_nullstream< char, char_traits< char >> nullstream
Definition: ROL_Stream.hpp:38
int main(int argc, char *argv[])
constexpr auto dim