1 // @HEADER
2 // ************************************************************************
3 //
4 // Rapid Optimization Library (ROL) Package
5 // Copyright (2014) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact lead developers:
38 // Drew Kouri (dpkouri@sandia.gov) and
39 // Denis Ridzal (dridzal@sandia.gov)
40 //
41 // ************************************************************************
42 // @HEADER
43
44 /*! \file example_03.cpp
45 \brief Shows how to solve the stochastic Navier-Stokes problem.
46 */
47
48 #include "Teuchos_Comm.hpp"
49 #ifdef HAVE_MPI
50 #include "Teuchos_DefaultMpiComm.hpp"
51 #endif
52 #include "ROL_Stream.hpp"
53 #include "Teuchos_GlobalMPISession.hpp"
54 #include "Teuchos_XMLParameterListHelpers.hpp"
55
56 #include "Tpetra_Core.hpp"
57 #include "Tpetra_Version.hpp"
58
59 #include <iostream>
60 #include <algorithm>
61 //#include <fenv.h>
62
63 #include "ROL_Bounds.hpp"
64 #include "ROL_Reduced_Objective_SimOpt.hpp"
65 #include "ROL_MonteCarloGenerator.hpp"
66 #include "ROL_OptimizationSolver.hpp"
67 #include "ROL_TpetraTeuchosBatchManager.hpp"
68
69 #include "../TOOLS/meshmanager.hpp"
70 #include "../TOOLS/pdeconstraint.hpp"
71 #include "../TOOLS/pdeobjective.hpp"
72 #include "../TOOLS/pdevector.hpp"
73 #include "../TOOLS/batchmanager.hpp"
74 #include "pde_navier-stokes.hpp"
75 #include "obj_navier-stokes.hpp"
76
77 typedef double RealT;
78
79 template<class Real>
random(const Teuchos::Comm<int> & comm,const Real a=-1,const Real b=1)80 Real random(const Teuchos::Comm<int> &comm,
81 const Real a = -1, const Real b = 1) {
82 Real val(0), u(0);
83 if ( Teuchos::rank<int>(comm)==0 ) {
84 u = static_cast<Real>(rand())/static_cast<Real>(RAND_MAX);
85 val = (b-a)*u + a;
86 }
87 Teuchos::broadcast<int,Real>(comm,0,1,&val);
88 return val;
89 }
90
91 template<class Real>
setUpAndSolve(ROL::OptimizationProblem<Real> & opt,Teuchos::ParameterList & parlist,std::ostream & outStream)92 void setUpAndSolve(ROL::OptimizationProblem<Real> &opt,
93 Teuchos::ParameterList &parlist,
94 std::ostream &outStream) {
95 parlist.sublist("Step").set("Type","Trust Region");
96 ROL::OptimizationSolver<Real> solver(opt,parlist);
97 Teuchos::Time timer("Optimization Time", true);
98 solver.solve(outStream);
99 timer.stop();
100 outStream << "Total optimization time = " << timer.totalElapsedTime() << " seconds." << std::endl;
101 }
102
103 template<class Real>
print(ROL::Objective<Real> & obj,const ROL::Vector<Real> & z,ROL::SampleGenerator<Real> & sampler,const int ngsamp,const ROL::Ptr<const Teuchos::Comm<int>> & comm,const std::string & filename)104 void print(ROL::Objective<Real> &obj,
105 const ROL::Vector<Real> &z,
106 ROL::SampleGenerator<Real> &sampler,
107 const int ngsamp,
108 const ROL::Ptr<const Teuchos::Comm<int> > &comm,
109 const std::string &filename) {
110 Real tol(1e-8);
111 // Build objective function distribution
112 int nsamp = sampler.numMySamples();
113 std::vector<Real> myvalues(nsamp), myzerovec(nsamp, 0);
114 std::vector<double> gvalues(ngsamp), gzerovec(ngsamp, 0);
115 std::vector<Real> sample = sampler.getMyPoint(0);
116 int sdim = sample.size();
117 std::vector<std::vector<Real> > mysamples(sdim, myzerovec);
118 std::vector<std::vector<double> > gsamples(sdim, gzerovec);
119 for (int i = 0; i < nsamp; ++i) {
120 sample = sampler.getMyPoint(i);
121 obj.setParameter(sample);
122 myvalues[i] = static_cast<double>(obj.value(z,tol));
123 for (int j = 0; j < sdim; ++j) {
124 mysamples[j][i] = static_cast<double>(sample[j]);
125 }
126 }
127
128 // Send data to root processor
129 #ifdef HAVE_MPI
130 ROL::Ptr<const Teuchos::MpiComm<int> > mpicomm
131 = ROL::dynamicPtrCast<const Teuchos::MpiComm<int> >(comm);
132 int nproc = Teuchos::size<int>(*mpicomm);
133 std::vector<int> sampleCounts(nproc, 0), sampleDispls(nproc, 0);
134 MPI_Gather(&nsamp,1,MPI_INT,&sampleCounts[0],1,MPI_INT,0,*(mpicomm->getRawMpiComm())());
135 for (int i = 1; i < nproc; ++i) {
136 sampleDispls[i] = sampleDispls[i-1] + sampleCounts[i-1];
137 }
138 MPI_Gatherv(&myvalues[0],nsamp,MPI_DOUBLE,&gvalues[0],&sampleCounts[0],&sampleDispls[0],MPI_DOUBLE,0,*(mpicomm->getRawMpiComm())());
139 for (int j = 0; j < sdim; ++j) {
140 MPI_Gatherv(&mysamples[j][0],nsamp,MPI_DOUBLE,&gsamples[j][0],&sampleCounts[0],&sampleDispls[0],MPI_DOUBLE,0,*(mpicomm->getRawMpiComm())());
141 }
142 #else
143 gvalues.assign(myvalues.begin(),myvalues.end());
144 for (int j = 0; j < sdim; ++j) {
145 gsamples[j].assign(mysamples[j].begin(),mysamples[j].end());
146 }
147 #endif
148
149 // Print
150 int rank = Teuchos::rank<int>(*comm);
151 if ( rank==0 ) {
152 std::ofstream file;
153 file.open(filename);
154 file << std::scientific << std::setprecision(15);
155 for (int i = 0; i < ngsamp; ++i) {
156 for (int j = 0; j < sdim; ++j) {
157 file << std::setw(25) << std::left << gsamples[j][i];
158 }
159 file << std::setw(25) << std::left << gvalues[i] << std::endl;
160 }
161 file.close();
162 }
163 }
164
main(int argc,char * argv[])165 int main(int argc, char *argv[]) {
166 // feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
167
168 // This little trick lets us print to std::cout only if a (dummy) command-line argument is provided.
169 int iprint = argc - 1;
170 ROL::Ptr<std::ostream> outStream;
171 ROL::nullstream bhs; // outputs nothing
172
173 /*** Initialize communicator. ***/
174 Teuchos::GlobalMPISession mpiSession (&argc, &argv, &bhs);
175 ROL::Ptr<const Teuchos::Comm<int> > comm
176 = Tpetra::getDefaultComm();
177 ROL::Ptr<const Teuchos::Comm<int> > serial_comm
178 = ROL::makePtr<Teuchos::SerialComm<int>>();
179 const int myRank = comm->getRank();
180 if ((iprint > 0) && (myRank == 0)) {
181 outStream = ROL::makePtrFromRef(std::cout);
182 }
183 else {
184 outStream = ROL::makePtrFromRef(bhs);
185 }
186 int errorFlag = 0;
187
188 // *** Example body.
189 try {
190
191 /*** Read in XML input ***/
192 std::string filename = "input_ex03.xml";
193 Teuchos::RCP<Teuchos::ParameterList> parlist = Teuchos::rcp( new Teuchos::ParameterList() );
194 Teuchos::updateParametersFromXmlFile( filename, parlist.ptr() );
195
196 // Problem dimensions
197 const int stochDim = 2;
198 const RealT one(1);
199
200 /*************************************************************************/
201 /***************** BUILD GOVERNING PDE ***********************************/
202 /*************************************************************************/
203 /*** Initialize main data structure. ***/
204 ROL::Ptr<MeshManager<RealT> > meshMgr
205 = ROL::makePtr<MeshManager_BackwardFacingStepChannel<RealT>>(*parlist);
206 // Initialize PDE describing advection-diffusion equation
207 ROL::Ptr<PDE_NavierStokes<RealT> > pde
208 = ROL::makePtr<PDE_NavierStokes<RealT>>(*parlist);
209 ROL::Ptr<ROL::Constraint_SimOpt<RealT> > con
210 = ROL::makePtr<PDE_Constraint<RealT>>(pde,meshMgr,serial_comm,*parlist,*outStream);
211 // Cast the constraint and get the assembler.
212 ROL::Ptr<PDE_Constraint<RealT> > pdecon
213 = ROL::dynamicPtrCast<PDE_Constraint<RealT> >(con);
214 ROL::Ptr<Assembler<RealT> > assembler = pdecon->getAssembler();
215 assembler->printMeshData(*outStream);
216 con->setSolveParameters(*parlist);
217
218 /*************************************************************************/
219 /***************** BUILD VECTORS *****************************************/
220 /*************************************************************************/
221 ROL::Ptr<Tpetra::MultiVector<> > u_ptr = assembler->createStateVector();
222 ROL::Ptr<Tpetra::MultiVector<> > p_ptr = assembler->createStateVector();
223 ROL::Ptr<Tpetra::MultiVector<> > du_ptr = assembler->createStateVector();
224 u_ptr->randomize(); //u_ptr->putScalar(static_cast<RealT>(1));
225 p_ptr->randomize(); //p_ptr->putScalar(static_cast<RealT>(1));
226 du_ptr->randomize(); //du_ptr->putScalar(static_cast<RealT>(0));
227 ROL::Ptr<ROL::Vector<RealT> > up
228 = ROL::makePtr<PDE_PrimalSimVector<RealT>>(u_ptr,pde,assembler,*parlist);
229 ROL::Ptr<ROL::Vector<RealT> > pp
230 = ROL::makePtr<PDE_PrimalSimVector<RealT>>(p_ptr,pde,assembler,*parlist);
231 ROL::Ptr<ROL::Vector<RealT> > dup
232 = ROL::makePtr<PDE_PrimalSimVector<RealT>>(du_ptr,pde,assembler,*parlist);
233 // Create residual vectors
234 ROL::Ptr<Tpetra::MultiVector<> > r_ptr = assembler->createResidualVector();
235 r_ptr->randomize(); //r_ptr->putScalar(static_cast<RealT>(1));
236 ROL::Ptr<ROL::Vector<RealT> > rp
237 = ROL::makePtr<PDE_DualSimVector<RealT>>(r_ptr,pde,assembler,*parlist);
238 // Create control vector and set to ones
239 ROL::Ptr<Tpetra::MultiVector<> > z_ptr = assembler->createControlVector();
240 ROL::Ptr<Tpetra::MultiVector<> > dz_ptr = assembler->createControlVector();
241 ROL::Ptr<Tpetra::MultiVector<> > yz_ptr = assembler->createControlVector();
242 z_ptr->randomize(); z_ptr->putScalar(static_cast<RealT>(0));
243 dz_ptr->randomize(); //dz_ptr->putScalar(static_cast<RealT>(0));
244 yz_ptr->randomize(); //yz_ptr->putScalar(static_cast<RealT>(0));
245 ROL::Ptr<ROL::TpetraMultiVector<RealT> > zpde
246 = ROL::makePtr<PDE_PrimalOptVector<RealT>>(z_ptr,pde,assembler,*parlist);
247 ROL::Ptr<ROL::TpetraMultiVector<RealT> > dzpde
248 = ROL::makePtr<PDE_PrimalOptVector<RealT>>(dz_ptr,pde,assembler,*parlist);
249 ROL::Ptr<ROL::TpetraMultiVector<RealT> > yzpde
250 = ROL::makePtr<PDE_PrimalOptVector<RealT>>(yz_ptr,pde,assembler,*parlist);
251 ROL::Ptr<ROL::Vector<RealT> > zp = ROL::makePtr<PDE_OptVector<RealT>>(zpde);
252 ROL::Ptr<ROL::Vector<RealT> > dzp = ROL::makePtr<PDE_OptVector<RealT>>(dzpde);
253 ROL::Ptr<ROL::Vector<RealT> > yzp = ROL::makePtr<PDE_OptVector<RealT>>(yzpde);
254 // Create ROL SimOpt vectors
255 ROL::Vector_SimOpt<RealT> x(up,zp);
256 ROL::Vector_SimOpt<RealT> d(dup,dzp);
257
258 /*************************************************************************/
259 /***************** BUILD COST FUNCTIONAL *********************************/
260 /*************************************************************************/
261 std::vector<ROL::Ptr<QoI<RealT> > > qoi_vec(2,ROL::nullPtr);
262 qoi_vec[0] = ROL::makePtr<QoI_State_NavierStokes<RealT>>(*parlist,
263 pde->getVelocityFE(),
264 pde->getPressureFE(),
265 pde->getFieldHelper());
266 qoi_vec[1] = ROL::makePtr<QoI_L2Penalty_NavierStokes<RealT>>(pde->getVelocityFE(),
267 pde->getPressureFE(),
268 pde->getVelocityBdryFE(),
269 pde->getBdryCellLocIds(),
270 pde->getFieldHelper());
271 ROL::Ptr<StdObjective_NavierStokes<RealT> > std_obj
272 = ROL::makePtr<StdObjective_NavierStokes<RealT>>(*parlist);
273 ROL::Ptr<ROL::Objective_SimOpt<RealT> > obj
274 = ROL::makePtr<PDE_Objective<RealT>>(qoi_vec,std_obj,assembler);
275 ROL::Ptr<ROL::Objective_SimOpt<RealT> > objState
276 = ROL::makePtr<PDE_Objective<RealT>>(qoi_vec[0],assembler);
277 ROL::Ptr<ROL::Objective_SimOpt<RealT> > objCtrl
278 = ROL::makePtr<PDE_Objective<RealT>>(qoi_vec[1],assembler);
279 ROL::Ptr<ROL::Reduced_Objective_SimOpt<RealT> > objRed
280 = ROL::makePtr<ROL::Reduced_Objective_SimOpt<RealT>>(obj, con, up, zp, pp, true, false);
281
282 /*************************************************************************/
283 /***************** BUILD BOUND CONSTRAINT ********************************/
284 /*************************************************************************/
285 ROL::Ptr<Tpetra::MultiVector<> > zlo_ptr = assembler->createControlVector();
286 ROL::Ptr<Tpetra::MultiVector<> > zhi_ptr = assembler->createControlVector();
287 zlo_ptr->putScalar(static_cast<RealT>(0));
288 zhi_ptr->putScalar(ROL::ROL_INF<RealT>());
289 ROL::Ptr<ROL::TpetraMultiVector<RealT> > zlopde
290 = ROL::makePtr<PDE_PrimalOptVector<RealT>>(zlo_ptr,pde,assembler,*parlist);
291 ROL::Ptr<ROL::TpetraMultiVector<RealT> > zhipde
292 = ROL::makePtr<PDE_PrimalOptVector<RealT>>(zhi_ptr,pde,assembler,*parlist);
293 ROL::Ptr<ROL::Vector<RealT> > zlop = ROL::makePtr<PDE_OptVector<RealT>>(zlopde);
294 ROL::Ptr<ROL::Vector<RealT> > zhip = ROL::makePtr<PDE_OptVector<RealT>>(zhipde);
295 ROL::Ptr<ROL::BoundConstraint<RealT> > bnd
296 = ROL::makePtr<ROL::Bounds<RealT>>(zlop,zhip);
297 bool useBounds = parlist->sublist("Problem").get("Use bounds", false);
298 if (!useBounds) {
299 bnd->deactivate();
300 }
301
302 /*************************************************************************/
303 /***************** BUILD SAMPLER *****************************************/
304 /*************************************************************************/
305 int nsamp = parlist->sublist("Problem").get("Number of samples",100);
306 int nsamp_dist = parlist->sublist("Problem").get("Number of output samples",100);
307 std::vector<RealT> tmp = {-one,one};
308 std::vector<std::vector<RealT> > bounds(stochDim,tmp);
309 ROL::Ptr<ROL::BatchManager<RealT> > bman
310 = ROL::makePtr<PDE_OptVector_BatchManager<RealT>>(comm);
311 ROL::Ptr<ROL::SampleGenerator<RealT> > sampler
312 = ROL::makePtr<ROL::MonteCarloGenerator<RealT>>(nsamp,bounds,bman);
313 ROL::Ptr<ROL::SampleGenerator<RealT> > sampler_dist
314 = ROL::makePtr<ROL::MonteCarloGenerator<RealT>>(nsamp_dist,bounds,bman);
315
316 /*************************************************************************/
317 /***************** BUILD STOCHASTIC PROBLEM ******************************/
318 /*************************************************************************/
319 ROL::Ptr<ROL::OptimizationProblem<RealT> > opt;
320 std::vector<RealT> ctrl;
321 std::vector<RealT> var;
322
323 Teuchos::Array<RealT> alphaArray
324 = Teuchos::getArrayFromStringParameter<RealT>(parlist->sublist("Problem"),"Confidence Levels");
325 std::vector<RealT> alpha = alphaArray.toVector();
326 std::sort(alpha.begin(),alpha.end());
327 int N = alpha.size();
328
329 /*************************************************************************/
330 /***************** SOLVE RISK NEUTRAL ************************************/
331 /*************************************************************************/
332 RealT tol(1e-8);
333 bool alphaZero = (alpha[0] == static_cast<RealT>(0));
334 if ( alphaZero ) {
335 alpha.erase(alpha.begin()); --N;
336 // Solve.
337 parlist->sublist("SOL").set("Stochastic Component Type","Risk Neutral");
338 opt = ROL::makePtr<ROL::OptimizationProblem<RealT>>(objRed,zp,bnd);
339 parlist->sublist("SOL").set("Initial Statistic",one);
340 opt->setStochasticObjective(*parlist,sampler);
341 setUpAndSolve<RealT>(*opt,*parlist,*outStream);
342 // Output.
343 ctrl.push_back(objCtrl->value(*up,*zp,tol));
344 var.push_back(opt->getSolutionStatistic());
345 std::string CtrlRN = "control_RN.txt";
346 pdecon->outputTpetraVector(z_ptr,CtrlRN);
347 std::string ObjRN = "obj_samples_RN.txt";
348 print<RealT>(*objRed,*zp,*sampler_dist,nsamp_dist,comm,ObjRN);
349 }
350
351 /*************************************************************************/
352 /***************** SOLVE MEAN PLUS CVAR **********************************/
353 /*************************************************************************/
354 parlist->sublist("SOL").set("Stochastic Component Type","Risk Averse");
355 parlist->sublist("SOL").sublist("Risk Measure").set("Name","CVaR");
356 parlist->sublist("SOL").sublist("Risk Measure").sublist("CVaR").set("Convex Combination Parameter",0.0);
357 parlist->sublist("SOL").sublist("Risk Measure").sublist("CVaR").set("Smoothing Parameter",1e-4);
358 parlist->sublist("SOL").sublist("Risk Measure").sublist("CVaR").sublist("Distribution").set("Name","Parabolic");
359 parlist->sublist("SOL").sublist("Risk Measure").sublist("CVaR").sublist("Distribution").sublist("Parabolic").set("Lower Bound",0.0);
360 parlist->sublist("SOL").sublist("Risk Measure").sublist("CVaR").sublist("Distribution").sublist("Parabolic").set("Upper Bound",1.0);
361 for (int i = 0; i < N; ++i) {
362 // Solve.
363 parlist->sublist("SOL").sublist("Risk Measure").sublist("CVaR").set("Confidence Level",alpha[i]);
364 opt = ROL::makePtr<ROL::OptimizationProblem<RealT>>(objRed,zp,bnd);
365 parlist->sublist("SOL").set("Initial Statistic",one);
366 opt->setStochasticObjective(*parlist,sampler);
367 setUpAndSolve<RealT>(*opt,*parlist,*outStream);
368 // Output.
369 ctrl.push_back(objCtrl->value(*up,*zp,tol));
370 var.push_back(opt->getSolutionStatistic());
371 std::stringstream nameCtrl;
372 nameCtrl << "control_CVaR_" << i+1 << ".txt";
373 pdecon->outputTpetraVector(z_ptr,nameCtrl.str().c_str());
374 std::stringstream nameObj;
375 nameObj << "obj_samples_CVaR_" << i+1 << ".txt";
376 print<RealT>(*objRed,*zp,*sampler_dist,nsamp_dist,comm,nameObj.str());
377 }
378
379 /*************************************************************************/
380 /***************** PRINT CONTROL OBJ AND VAR *****************************/
381 /*************************************************************************/
382 const int rank = Teuchos::rank<int>(*comm);
383 if ( rank==0 ) {
384 std::stringstream nameCTRL, nameVAR;
385 nameCTRL << "ctrl.txt";
386 nameVAR << "var.txt";
387 std::ofstream fileCTRL, fileVAR;
388 fileCTRL.open(nameCTRL.str());
389 fileVAR.open(nameVAR.str());
390 fileCTRL << std::scientific << std::setprecision(15);
391 fileVAR << std::scientific << std::setprecision(15);
392 int size = var.size();
393 for (int i = 0; i < size; ++i) {
394 fileCTRL << std::setw(25) << std::left << ctrl[i] << std::endl;
395 fileVAR << std::setw(25) << std::left << var[i] << std::endl;
396 }
397 fileCTRL.close();
398 fileVAR.close();
399 }
400
401 // Get a summary from the time monitor.
402 Teuchos::TimeMonitor::summarize();
403 }
404 catch (std::logic_error& err) {
405 *outStream << err.what() << "\n";
406 errorFlag = -1000;
407 }; // end try
408
409 if (errorFlag != 0)
410 std::cout << "End Result: TEST FAILED\n";
411 else
412 std::cout << "End Result: TEST PASSED\n";
413
414 return 0;
415 }
416