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 Poisson-Boltzmann problem.
46 */
47
48 #include "Teuchos_Comm.hpp"
49 #include "Teuchos_GlobalMPISession.hpp"
50 #include "Teuchos_XMLParameterListHelpers.hpp"
51
52 #include "Tpetra_Core.hpp"
53 #include "Tpetra_Version.hpp"
54
55 #include <iostream>
56 #include <algorithm>
57
58 #include "ROL_Stream.hpp"
59 #include "ROL_OptimizationSolver.hpp"
60 #include "ROL_Reduced_Objective_SimOpt.hpp"
61 #include "ROL_BoundConstraint_SimOpt.hpp"
62 #include "ROL_Bounds.hpp"
63
64 #include "../TOOLS/meshmanager.hpp"
65 #include "../TOOLS/pdeconstraint.hpp"
66 #include "../TOOLS/pdeobjective.hpp"
67 #include "../TOOLS/pdevector.hpp"
68 #include "pde_allen_cahn.hpp"
69 #include "obj_allen_cahn.hpp"
70
71 typedef double RealT;
72
main(int argc,char * argv[])73 int main(int argc, char *argv[]) {
74 // This little trick lets us print to std::cout only if a (dummy) command-line argument is provided.
75 int iprint = argc - 1;
76 ROL::Ptr<std::ostream> outStream;
77 ROL::nullstream bhs; // outputs nothing
78
79 /*** Initialize communicator. ***/
80 Teuchos::GlobalMPISession mpiSession (&argc, &argv, &bhs);
81 ROL::Ptr<const Teuchos::Comm<int> > comm
82 = Tpetra::getDefaultComm();
83 const int myRank = comm->getRank();
84 if ((iprint > 0) && (myRank == 0)) {
85 outStream = ROL::makePtrFromRef(std::cout);
86 }
87 else {
88 outStream = ROL::makePtrFromRef(bhs);
89 }
90 int errorFlag = 0;
91
92 // *** Example body.
93 try {
94
95 /*** Read in XML input ***/
96 std::string filename = "input.xml";
97 Teuchos::RCP<Teuchos::ParameterList> parlist = Teuchos::rcp( new Teuchos::ParameterList() );
98 Teuchos::updateParametersFromXmlFile( filename, parlist.ptr() );
99
100 /*** Initialize main data structure. ***/
101 ROL::Ptr<MeshManager<RealT> > meshMgr
102 = ROL::makePtr<MeshManager_Rectangle<RealT>>(*parlist);
103 // Initialize PDE describe Poisson's equation
104 ROL::Ptr<PDE_Allen_Cahn<RealT> > pde
105 = ROL::makePtr<PDE_Allen_Cahn<RealT>>(*parlist);
106 ROL::Ptr<PDE_Constraint<RealT> > con
107 = ROL::makePtr<PDE_Constraint<RealT>>(pde,meshMgr,comm,*parlist,*outStream);
108 const ROL::Ptr<Assembler<RealT> > assembler = con->getAssembler();
109 assembler->printMeshData(*outStream);
110 con->setSolveParameters(*parlist);
111
112 /*************************************************************************/
113 /***************** BUILD VECTORS *****************************************/
114 /*************************************************************************/
115 ROL::Ptr<Tpetra::MultiVector<> > u_ptr = assembler->createStateVector();
116 ROL::Ptr<Tpetra::MultiVector<> > p_ptr = assembler->createStateVector();
117 ROL::Ptr<Tpetra::MultiVector<> > r_ptr = assembler->createResidualVector();
118 ROL::Ptr<Tpetra::MultiVector<> > z_ptr = assembler->createControlVector();
119 ROL::Ptr<ROL::Vector<RealT> > up, pp, rp, zp;
120 u_ptr->randomize(); //u_ptr->putScalar(static_cast<RealT>(1));
121 p_ptr->randomize(); //p_ptr->putScalar(static_cast<RealT>(1));
122 r_ptr->randomize(); //r_ptr->putScalar(static_cast<RealT>(1));
123 z_ptr->randomize(); //z_ptr->putScalar(static_cast<RealT>(1));
124 up = ROL::makePtr<PDE_PrimalSimVector<RealT>>(u_ptr,pde,assembler,*parlist);
125 pp = ROL::makePtr<PDE_PrimalSimVector<RealT>>(p_ptr,pde,assembler,*parlist);
126 rp = ROL::makePtr<PDE_DualSimVector<RealT>>(r_ptr,pde,assembler,*parlist);
127 zp = ROL::makePtr<PDE_PrimalOptVector<RealT>>(z_ptr,pde,assembler,*parlist);
128
129 /*************************************************************************/
130 /***************** BUILD COST FUNCTIONAL *********************************/
131 /*************************************************************************/
132 // Initialize quadratic objective function
133 std::vector<ROL::Ptr<QoI<RealT> > > qoi_vec(2,ROL::nullPtr);
134 qoi_vec[0] = ROL::makePtr<QoI_State_Cost_Allen_Cahn<RealT>>(pde->getFE());
135 qoi_vec[1] = ROL::makePtr<QoI_Control_Cost_Allen_Cahn<RealT>>(pde->getFE(), pde->getBdryFE(), pde->getBdryCellLocIds());
136 std::vector<RealT> weights(2);
137 weights[0] = static_cast<RealT>(1);
138 weights[1] = parlist->sublist("Problem").get("Control Penalty Parameter", 1e-4);
139 // Build full-space objective
140 ROL::Ptr<PDE_Objective<RealT> > obj
141 = ROL::makePtr<PDE_Objective<RealT>>(qoi_vec,weights,assembler);
142 // Build reduced-space objective
143 bool storage = parlist->sublist("Problem").get("Use state storage",true);
144 ROL::Ptr<ROL::SimController<RealT> > stateStore
145 = ROL::makePtr<ROL::SimController<RealT>>();
146 ROL::Ptr<ROL::Reduced_Objective_SimOpt<RealT> > robj
147 = ROL::makePtr<ROL::Reduced_Objective_SimOpt<RealT>>(
148 obj,con,stateStore,up,zp,pp,storage);
149
150 /*************************************************************************/
151 /***************** BUILD BOUND CONSTRAINT ********************************/
152 /*************************************************************************/
153 // Control bounds
154 ROL::Ptr<Tpetra::MultiVector<> > zlo_ptr = assembler->createControlVector();
155 ROL::Ptr<Tpetra::MultiVector<> > zhi_ptr = assembler->createControlVector();
156 RealT lo = parlist->sublist("Problem").get("Lower Bound",0.0);
157 RealT hi = parlist->sublist("Problem").get("Upper Bound",0.0);
158 zlo_ptr->putScalar(lo); zhi_ptr->putScalar(hi);
159 ROL::Ptr<ROL::Vector<RealT> > zlop, zhip;
160 zlop = ROL::makePtr<PDE_PrimalOptVector<RealT>>(zlo_ptr,pde,assembler);
161 zhip = ROL::makePtr<PDE_PrimalOptVector<RealT>>(zhi_ptr,pde,assembler);
162 ROL::Ptr<ROL::BoundConstraint<RealT> > zbnd
163 = ROL::makePtr<ROL::Bounds<RealT>>(zlop,zhip);
164 bool deactivate = parlist->sublist("Problem").get("Deactivate Bound Constraints",false);
165 if (deactivate) {
166 zbnd->deactivate();
167 }
168 // State bounds
169 ROL::Ptr<Tpetra::MultiVector<> > ulo_ptr = assembler->createStateVector();
170 ROL::Ptr<Tpetra::MultiVector<> > uhi_ptr = assembler->createStateVector();
171 ulo_ptr->putScalar(ROL::ROL_NINF<RealT>()); uhi_ptr->putScalar(ROL::ROL_INF<RealT>());
172 ROL::Ptr<ROL::Vector<RealT> > ulop
173 = ROL::makePtr<PDE_PrimalSimVector<RealT>>(ulo_ptr,pde,assembler);
174 ROL::Ptr<ROL::Vector<RealT> > uhip
175 = ROL::makePtr<PDE_PrimalSimVector<RealT>>(uhi_ptr,pde,assembler);
176 ROL::Ptr<ROL::BoundConstraint<RealT> > ubnd
177 = ROL::makePtr<ROL::Bounds<RealT>>(ulop,uhip);
178 ubnd->deactivate();
179
180 // SimOpt bounds
181 ROL::Ptr<ROL::BoundConstraint<RealT> > bnd
182 = ROL::makePtr<ROL::BoundConstraint_SimOpt<RealT>>(ubnd,zbnd);
183
184 // Create ROL SimOpt vectors
185 ROL::Ptr<Tpetra::MultiVector<> > du_ptr = assembler->createStateVector();
186 ROL::Ptr<Tpetra::MultiVector<> > dz_ptr = assembler->createControlVector();
187 du_ptr->randomize(); //du_ptr->putScalar(static_cast<RealT>(0));
188 dz_ptr->randomize(); //dz_ptr->putScalar(static_cast<RealT>(1));
189 ROL::Ptr<ROL::Vector<RealT> > dup, dzp;
190 dup = ROL::makePtr<PDE_PrimalSimVector<RealT>>(du_ptr,pde,assembler,*parlist);
191 dzp = ROL::makePtr<PDE_PrimalOptVector<RealT>>(dz_ptr,pde,assembler,*parlist);
192
193 ROL::Vector_SimOpt<RealT> x(up,zp);
194 ROL::Vector_SimOpt<RealT> d(dup,dzp);
195
196 /*************************************************************************/
197 /***************** RUN VECTOR AND DERIVATIVE CHECKS **********************/
198 /*************************************************************************/
199 bool checkDeriv = parlist->sublist("Problem").get("Check Derivatives",false);
200 if ( checkDeriv ) {
201 *outStream << "\n\nCheck Gradient of Full Objective Function\n";
202 obj->checkGradient(x,d,true,*outStream);
203 *outStream << "\n\nCheck Hessian of Full Objective Function\n";
204 obj->checkHessVec(x,d,true,*outStream);
205 *outStream << "\n\nCheck Full Jacobian of PDE Constraint\n";
206 con->checkApplyJacobian(x,d,*rp,true,*outStream);
207 *outStream << "\n\nCheck Jacobian_1 of PDE Constraint\n";
208 con->checkApplyJacobian_1(*up,*zp,*dup,*rp,true,*outStream);
209 *outStream << "\n\nCheck Jacobian_2 of PDE Constraint\n";
210 con->checkApplyJacobian_2(*up,*zp,*dzp,*rp,true,*outStream);
211 *outStream << "\n\nCheck Full Hessian of PDE Constraint\n";
212 con->checkApplyAdjointHessian(x,*pp,d,x,true,*outStream);
213 *outStream << "\n";
214 con->checkAdjointConsistencyJacobian(*dup,d,x,true,*outStream);
215 *outStream << "\n";
216 con->checkInverseJacobian_1(*rp,*dup,*up,*zp,true,*outStream);
217 *outStream << "\n";
218 *outStream << "\n\nCheck Gradient of Reduced Objective Function\n";
219 robj->checkGradient(*zp,*dzp,true,*outStream);
220 *outStream << "\n\nCheck Hessian of Reduced Objective Function\n";
221 robj->checkHessVec(*zp,*dzp,true,*outStream);
222 }
223
224 bool useFullSpace = parlist->sublist("Problem").get("Full space",false);
225 ROL::Ptr<ROL::Algorithm<RealT> > algo;
226 if ( useFullSpace ) {
227 ROL::OptimizationProblem<RealT> optProb(obj, makePtrFromRef(x), bnd, con, pp);
228 ROL::OptimizationSolver<RealT> optSolver(optProb, *parlist);
229 optSolver.solve(*outStream);
230 }
231 else {
232 ROL::Ptr<ROL::Step<RealT>>
233 step = ROL::makePtr<ROL::TrustRegionStep<RealT>>(*parlist);
234 ROL::Ptr<ROL::StatusTest<RealT>>
235 status = ROL::makePtr<ROL::StatusTest<RealT>>(*parlist);
236 ROL::Algorithm<RealT> algo(step,status,false);
237 algo.run(*zp,*robj,*zbnd,true,*outStream);
238 }
239
240 // Output.
241 RealT tol(1.e-8);
242 con->solve(*rp,*up,*zp,tol);
243 con->outputTpetraVector(u_ptr,"state.txt");
244 con->outputTpetraVector(z_ptr,"control.txt");
245
246 Teuchos::Array<RealT> res(1,0);
247 con->value(*rp,*up,*zp,tol);
248 r_ptr->norm2(res.view(0,1));
249 *outStream << "Residual Norm: " << res[0] << std::endl;
250 errorFlag += (res[0] > 1.e-6 ? 1 : 0);
251 }
252 catch (std::logic_error& err) {
253 *outStream << err.what() << "\n";
254 errorFlag = -1000;
255 }; // end try
256
257 if (errorFlag != 0)
258 std::cout << "End Result: TEST FAILED\n";
259 else
260 std::cout << "End Result: TEST PASSED\n";
261
262 return 0;
263 }
264