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_01.cpp
45 \brief Shows how to solve the mother problem of PDE-constrained optimization:
46 */
47
48 #include "Teuchos_Comm.hpp"
49 #include "ROL_Stream.hpp"
50 #include "Teuchos_GlobalMPISession.hpp"
51 #include "Teuchos_XMLParameterListHelpers.hpp"
52
53 #include "Tpetra_Core.hpp"
54 #include "Tpetra_Version.hpp"
55
56 #include "ROL_Algorithm.hpp"
57 #include "ROL_AugmentedLagrangianStep.hpp"
58 //#include "ROL_MoreauYosidaPenaltyStep.hpp"
59 #include "ROL_TrustRegionStep.hpp"
60 #include "ROL_CompositeStep.hpp"
61 #include "ROL_Reduced_Objective_SimOpt.hpp"
62 #include "ROL_ScaledTpetraMultiVector.hpp"
63
64 #include "ROL_ScaledStdVector.hpp"
65
66 #include "ROL_Bounds.hpp"
67 #include "ROL_AugmentedLagrangian.hpp"
68 #include "ROL_Algorithm.hpp"
69
70 #include <iostream>
71 #include <algorithm>
72
73 #include "data.hpp"
74 #include "filter.hpp"
75 #include "objective.hpp"
76 #include "constraint.hpp"
77 #include "volume_constraint.hpp"
78
79 //#include <fenv.h>
80
81 typedef double RealT;
82
createTpetraVector(const ROL::Ptr<const Tpetra::Map<>> & map)83 ROL::Ptr<Tpetra::MultiVector<> > createTpetraVector(const ROL::Ptr<const Tpetra::Map<> > &map) {
84 return ROL::makePtr<Tpetra::MultiVector<>>(map, 1, true);
85 }
86
main(int argc,char * argv[])87 int main(int argc, char *argv[]) {
88 //feenableexcept(FE_ALL_EXCEPT & ~FE_INEXACT);
89
90 // This little trick lets us print to std::cout only if a (dummy) command-line argument is provided.
91 int iprint = argc - 1;
92 ROL::Ptr<std::ostream> outStream;
93 ROL::nullstream bhs; // outputs nothing
94
95 /*** Initialize communicator. ***/
96 Teuchos::GlobalMPISession mpiSession (&argc, &argv, &bhs);
97 ROL::Ptr<const Teuchos::Comm<int> > comm = Tpetra::getDefaultComm();
98 const int myRank = comm->getRank();
99 if ((iprint > 0) && (myRank == 0)) {
100 outStream = ROL::makePtrFromRef(std::cout);
101 }
102 else {
103 outStream = ROL::makePtrFromRef(bhs);
104 }
105
106 int errorFlag = 0;
107
108 // *** Example body.
109 try {
110
111 /*** Read in XML input ***/
112 std::string filename = "input.xml";
113 Teuchos::RCP<Teuchos::ParameterList> parlist = Teuchos::rcp( new Teuchos::ParameterList() );
114 Teuchos::updateParametersFromXmlFile( filename, parlist.ptr() );
115
116 /*** Initialize main data structure. ***/
117 ROL::Ptr<ElasticitySIMPOperators<RealT> > data
118 = ROL::makePtr<ElasticitySIMPOperators<RealT>>(comm, parlist, outStream);
119 /*** Initialize density filter. ***/
120 ROL::Ptr<DensityFilter<RealT> > filter
121 = ROL::makePtr<DensityFilter<RealT>>(comm, parlist, outStream);
122 /*** Build vectors and dress them up as ROL vectors. ***/
123 ROL::Ptr<const Tpetra::Map<> > vecmap_u = data->getDomainMapA();
124 ROL::Ptr<const Tpetra::Map<> > vecmap_z = data->getCellMap();
125 ROL::Ptr<Tpetra::MultiVector<> > u_ptr = createTpetraVector(vecmap_u);
126 ROL::Ptr<Tpetra::MultiVector<> > z_ptr = createTpetraVector(vecmap_z);
127 ROL::Ptr<Tpetra::MultiVector<> > du_ptr = createTpetraVector(vecmap_u);
128 ROL::Ptr<Tpetra::MultiVector<> > dw_ptr = createTpetraVector(vecmap_u);
129 ROL::Ptr<Tpetra::MultiVector<> > dz_ptr = createTpetraVector(vecmap_z);
130 ROL::Ptr<Tpetra::MultiVector<> > dz2_ptr = createTpetraVector(vecmap_z);
131 ROL::Ptr<std::vector<RealT> > vc_ptr = ROL::makePtr<std::vector<RealT>>(1, 0);
132 ROL::Ptr<std::vector<RealT> > vc_lam_ptr = ROL::makePtr<std::vector<RealT>>(1, 0);
133 ROL::Ptr<std::vector<RealT> > vscale_ptr = ROL::makePtr<std::vector<RealT>>(1, 0);
134 // Set all values to 1 in u, z.
135 u_ptr->putScalar(1.0);
136 // Set z to gray solution.
137 RealT volFrac = parlist->sublist("ElasticityTopoOpt").get<RealT>("Volume Fraction");
138 z_ptr->putScalar(volFrac);
139 // Set scaling vector for constraint
140 RealT W = parlist->sublist("Geometry").get<RealT>("Width");
141 RealT H = parlist->sublist("Geometry").get<RealT>("Height");
142 RealT one(1), two(2);
143 (*vscale_ptr)[0] = one/std::pow(W*H*(one-volFrac),two);
144 // Set Scaling vector for density
145 bool useZscale = parlist->sublist("Problem").get<bool>("Use Scaled Density Vectors");
146 RealT densityScaling = parlist->sublist("Problem").get<RealT>("Density Scaling");
147 ROL::Ptr<Tpetra::MultiVector<> > scaleVec = createTpetraVector(vecmap_z);
148 scaleVec->putScalar(densityScaling);
149 if ( !useZscale ) {
150 scaleVec->putScalar(one);
151 }
152 ROL::Ptr<const Tpetra::Vector<> > zscale_ptr = scaleVec->getVector(0);
153
154 //test
155 /*data->updateMaterialDensity (z_ptr);
156 ROL::Ptr<Tpetra::MultiVector<RealT> > rhs
157 = ROL::makePtr<Tpetra::MultiVector<>> (data->getVecF()->getMap(), 1, true);
158 data->ApplyMatAToVec(rhs, u_ptr);
159 data->outputTpetraVector(rhs, "KU0.txt");
160 data->ApplyInverseJacobian1ToVec(u_ptr, rhs, false);
161 data->outputTpetraVector(u_ptr, "KKU0.txt");
162
163 data->ApplyJacobian1ToVec(rhs, u_ptr);
164 data->outputTpetraVector(rhs, "KU1.txt");
165 data->ApplyInverseJacobian1ToVec(u_ptr, rhs, false);
166 data->outputTpetraVector(u_ptr, "KKU1.txt");
167 */
168 //u_ptr->putScalar(1.0);
169 //z_ptr->putScalar(1.0);
170 // Randomize d vectors.
171 du_ptr->randomize(); //du_ptr->scale(0);
172 dw_ptr->randomize();
173 dz_ptr->randomize(); //dz_ptr->scale(0);
174 dz2_ptr->randomize();
175 // Create ROL::TpetraMultiVectors.
176 ROL::Ptr<ROL::Vector<RealT> > up = ROL::makePtr<ROL::TpetraMultiVector<RealT>>(u_ptr);
177 ROL::Ptr<ROL::Vector<RealT> > dup = ROL::makePtr<ROL::TpetraMultiVector<RealT>>(du_ptr);
178 ROL::Ptr<ROL::Vector<RealT> > dwp = ROL::makePtr<ROL::TpetraMultiVector<RealT>>(dw_ptr);
179 ROL::Ptr<ROL::Vector<RealT> > zp
180 = ROL::makePtr<ROL::PrimalScaledTpetraMultiVector<RealT>>(z_ptr,zscale_ptr);
181 ROL::Ptr<ROL::Vector<RealT> > dzp
182 = ROL::makePtr<ROL::PrimalScaledTpetraMultiVector<RealT>>(dz_ptr,zscale_ptr);
183 ROL::Ptr<ROL::Vector<RealT> > dz2p
184 = ROL::makePtr<ROL::PrimalScaledTpetraMultiVector<RealT>>(dz2_ptr,zscale_ptr);
185 ROL::Ptr<ROL::Vector<RealT> > vcp
186 = ROL::makePtr<ROL::PrimalScaledStdVector<RealT>>(vc_ptr,vscale_ptr);
187 ROL::Ptr<ROL::Vector<RealT> > vc_lamp
188 = ROL::makePtr<ROL::DualScaledStdVector<RealT>>(vc_lam_ptr,vscale_ptr);
189 // Create ROL SimOpt vectors.
190 ROL::Vector_SimOpt<RealT> x(up,zp);
191 ROL::Vector_SimOpt<RealT> d(dup,dzp);
192 ROL::Vector_SimOpt<RealT> d2(dwp,dz2p);
193
194 /*** Build objective function, constraint and reduced objective function. ***/
195 ROL::Ptr<ROL::Objective_SimOpt<RealT> > obj
196 = ROL::makePtr<Objective_PDEOPT_ElasticitySIMP<RealT>>(data, filter, parlist);
197 ROL::Ptr<ROL::Constraint_SimOpt<RealT> > con
198 = ROL::makePtr<EqualityConstraint_PDEOPT_ElasticitySIMP<RealT>>(data, filter, parlist);
199 ROL::Ptr<ROL::Reduced_Objective_SimOpt<RealT> > objReduced
200 = ROL::makePtr<ROL::Reduced_Objective_SimOpt<RealT>>(obj, con, up, zp, dwp);
201 ROL::Ptr<ROL::Constraint<RealT> > volcon
202 = ROL::makePtr<EqualityConstraint_PDEOPT_ElasticitySIMP_Volume<RealT>>(data, parlist);
203
204 /*** Build bound constraint ***/
205 ROL::Ptr<Tpetra::MultiVector<> > lo_ptr = ROL::makePtr<Tpetra::MultiVector<>>(vecmap_z, 1, true);
206 ROL::Ptr<Tpetra::MultiVector<> > hi_ptr = ROL::makePtr<Tpetra::MultiVector<>>(vecmap_z, 1, true);
207 lo_ptr->putScalar(0.0); hi_ptr->putScalar(1.0);
208 ROL::Ptr<ROL::Vector<RealT> > lop
209 = ROL::makePtr<ROL::PrimalScaledTpetraMultiVector<RealT>>(lo_ptr, zscale_ptr);
210 ROL::Ptr<ROL::Vector<RealT> > hip
211 = ROL::makePtr<ROL::PrimalScaledTpetraMultiVector<RealT>>(hi_ptr, zscale_ptr);
212 ROL::Ptr<ROL::BoundConstraint<RealT> > bnd = ROL::makePtr<ROL::Bounds<RealT>>(lop,hip);
213
214 /*** Check functional interface. ***/
215 *outStream << "Checking Objective:" << "\n";
216 obj->checkGradient(x,d,true,*outStream);
217 obj->checkHessVec(x,d,true,*outStream);
218 obj->checkHessSym(x,d,d2,true,*outStream);
219 *outStream << "Checking Constraint:" << "\n";
220 con->checkAdjointConsistencyJacobian(*dup,d,x,true,*outStream);
221 con->checkAdjointConsistencyJacobian_1(*dwp, *dup, *up, *zp, true, *outStream);
222 con->checkAdjointConsistencyJacobian_2(*dwp, *dzp, *up, *zp, true, *outStream);
223 con->checkInverseJacobian_1(*up,*up,*up,*zp,true,*outStream);
224 con->checkInverseAdjointJacobian_1(*up,*up,*up,*zp,true,*outStream);
225 con->checkApplyJacobian(x,d,*up,true,*outStream);
226 con->checkApplyAdjointHessian(x,*dup,d,x,true,*outStream);
227 *outStream << "Checking Reduced Objective:" << "\n";
228 objReduced->checkGradient(*zp,*dzp,true,*outStream);
229 objReduced->checkHessVec(*zp,*dzp,true,*outStream);
230 *outStream << "Checking Volume Constraint:" << "\n";
231 volcon->checkAdjointConsistencyJacobian(*vcp,*dzp,*zp,true,*outStream);
232 volcon->checkApplyJacobian(*zp,*dzp,*vcp,true,*outStream);
233 volcon->checkApplyAdjointHessian(*zp,*vcp,*dzp,*zp,true,*outStream);
234
235 /*** Run optimization ***/
236 ROL::Ptr<ROL::Step<RealT>>
237 step = ROL::makePtr<ROL::AugmentedLagrangianStep<RealT>>(*parlist);
238 ROL::Ptr<ROL::StatusTest<RealT>>
239 status = ROL::makePtr<ROL::ConstraintStatusTest<RealT>>(*parlist);
240 ROL::Algorithm<RealT> algo(step,status,false);
241 ROL::AugmentedLagrangian<RealT> augLag(objReduced,volcon,*vc_lamp,1.0,*zp,*vcp,*parlist);
242 algo.run(*zp,*vc_lamp,augLag,*volcon,*bnd,true,*outStream);
243 //ROL::Ptr<ROL::Step<RealT>>
244 // step = ROL::makePtr<ROL::MoreauYosidaPenaltyStep<RealT>>(*parlist);
245 //ROL::Ptr<ROL::StatusTest<RealT>>
246 // status = ROL::makePtr<ROL::ConstraintStatusTest<RealT>>(*parlist);
247 //ROL::MoreauYosidaPenalty<RealT> MYpen(objReduced,bnd,*zp,*parlist);
248 //ROL::Algorithm<RealT> algo(step,status,false);
249 //algo.run(*zp,*vc_lamp,MYpen,*volcon,*bnd,true,*outStream);
250
251 // new filter, for testing
252 /*parlist->sublist("Density Filter").set("Enable", true);
253 ROL::Ptr<DensityFilter<RealT> > testfilter
254 = ROL::makePtr<DensityFilter<RealT>>(comm, parlist, outStream);
255 ROL::Ptr<Tpetra::MultiVector<> > z_filtered_ptr = ROL::makePtr<Tpetra::MultiVector<>>(*z_ptr, Teuchos::Copy);
256 testfilter->apply(z_filtered_ptr, z_ptr);
257 ROL::Ptr<Tpetra::MultiVector<> > cm_ptr = data->getCellAreas();
258 ROL::Ptr<Tpetra::MultiVector<> > icm_ptr = ROL::makePtr<Tpetra::MultiVector<>>(*cm_ptr, Teuchos::Copy);
259 ROL::Ptr<Tpetra::MultiVector<> > zf_scaled_ptr = ROL::makePtr<Tpetra::MultiVector<>>(*z_ptr, Teuchos::Copy);
260 icm_ptr->reciprocal(*cm_ptr);
261 zf_scaled_ptr->elementWiseMultiply(1.0, *(icm_ptr->getVector(0)), *z_filtered_ptr, 0.0);
262 data->outputTpetraVector(zf_scaled_ptr, "density_filtered_scaled.txt");
263 data->outputTpetraVector(z_filtered_ptr, "density_filtered.txt");*/
264
265 data->outputTpetraVector(z_ptr, "density.txt");
266 data->outputTpetraVector(u_ptr, "state.txt");
267 data->outputTpetraVector(zscale_ptr, "weights.txt");
268
269 // Get a summary from the time monitor.
270 Teuchos::TimeMonitor::summarize();
271 }
272 catch (std::logic_error& err) {
273 *outStream << err.what() << "\n";
274 errorFlag = -1000;
275 }; // end try
276
277 if (errorFlag != 0)
278 std::cout << "End Result: TEST FAILED\n";
279 else
280 std::cout << "End Result: TEST PASSED\n";
281
282 return 0;
283 }
284