ROL
function/test_04.cpp
Go to the documentation of this file.
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
48#include "ROL_Stream.hpp"
49
50#include "Teuchos_GlobalMPISession.hpp"
51#include "Teuchos_Comm.hpp"
52#include "Teuchos_DefaultComm.hpp"
53#include "Teuchos_CommHelpers.hpp"
54
55#include <iostream>
56#include <fstream>
57#include <algorithm>
58
59#include "test_04.hpp"
60
61typedef double RealT;
68
69int main(int argc, char *argv[]) {
70
71 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
72
73 // This little trick lets us print to std::cout only if a (dummy) command-line argument is provided.
74 int iprint = argc - 1;
75 bool print = (iprint>0);
76 ROL::Ptr<std::ostream> outStream;
77 ROL::nullstream bhs; // outputs nothing
78 if (print)
79 outStream = ROL::makePtrFromRef(std::cout);
80 else
81 outStream = ROL::makePtrFromRef(bhs);
82
83 int errorFlag = 0;
84
85 // *** Example body.
86
87 try {
88 /*************************************************************************/
89 /************* INITIALIZE BURGERS FEM CLASS ******************************/
90 /*************************************************************************/
91 int nx = 512; // Set spatial discretization.
92 RealT nl = 1.0; // Nonlinearity parameter (1 = Burgers, 0 = linear).
93 RealT cH1 = 1.0; // Scale for derivative term in H1 norm.
94 RealT cL2 = 0.0; // Scale for mass term in H1 norm.
95 ROL::Ptr<BurgersFEM<RealT> > fem
96 = ROL::makePtr<BurgersFEM<RealT>>(nx,nl,cH1,cL2);
97 fem->test_inverse_mass(*outStream);
98 fem->test_inverse_H1(*outStream);
99 /*************************************************************************/
100 /************* INITIALIZE SIMOPT CONSTRAINT ******************************/
101 /*************************************************************************/
102 bool hess = true;
103 ROL::Ptr<ROL::Constraint_SimOpt<RealT> > con
104 = ROL::makePtr<Constraint_BurgersControl<RealT>>(fem,hess);
105 /*************************************************************************/
106 /************* INITIALIZE VECTOR STORAGE *********************************/
107 /*************************************************************************/
108 // INITIALIZE CONTROL VECTORS
109 ROL::Ptr<std::vector<RealT> > z_ptr
110 = ROL::makePtr<std::vector<RealT>>(nx+2, 0.0);
111 ROL::Ptr<ROL::Vector<RealT> > zp
112 = ROL::makePtr<PrimalControlVector>(z_ptr,fem);
113 // INITIALIZE STATE VECTORS
114 ROL::Ptr<std::vector<RealT> > u_ptr
115 = ROL::makePtr<std::vector<RealT>>(nx, 1.0);
116 ROL::Ptr<ROL::Vector<RealT> > up
117 = ROL::makePtr<PrimalStateVector>(u_ptr,fem);
118 // INITIALIZE CONSTRAINT VECTORS
119 ROL::Ptr<std::vector<RealT> > c_ptr
120 = ROL::makePtr<std::vector<RealT>>(nx, 1.0);
121 ROL::Ptr<ROL::Vector<RealT> > cp
122 = ROL::makePtr<PrimalConstraintVector>(c_ptr,fem);
123 /*************************************************************************/
124 /************* CHECK DERIVATIVES AND CONSISTENCY *************************/
125 /*************************************************************************/
126 RealT rnorm(0), cnorm(0);
127 ROL::ParameterList list;
128 RealT tol = std::sqrt(ROL::ROL_EPSILON<RealT>());
129 list.sublist("SimOpt").sublist("Solve").set("Output Iteration History",print);
130 list.sublist("SimOpt").sublist("Solve").set("Step Tolerance",ROL::ROL_EPSILON<RealT>());
131 // Newton
132 list.sublist("SimOpt").sublist("Solve").set("Absolute Residual Tolerance",tol);
133 list.sublist("SimOpt").sublist("Solve").set("Solver Type",0);
134 con->setSolveParameters(list);
135 u_ptr->assign(nx, 1.0);
136 con->solve(*cp,*up,*zp,tol);
137 rnorm = cp->norm();
138 con->value(*cp,*up,*zp,tol);
139 cnorm = cp->norm();
140 errorFlag += ((cnorm > tol) ? 1 : 0) + ((rnorm > tol) ? 1 : 0);
141 *outStream << std::scientific << std::setprecision(8);
142 *outStream << std::endl;
143 *outStream << "Test SimOpt solve at feasible (u,z):" << std::endl;
144 *outStream << " Solver Residual = " << rnorm << std::endl;
145 *outStream << " ||c(u,z)|| = " << cnorm;
146 *outStream << std::endl << std::endl;
147 // Levenberg-Marquardt
148 list.sublist("SimOpt").sublist("Solve").set("Absolute Residual Tolerance",1e-4*tol);
149 list.sublist("SimOpt").sublist("Solve").set("Solver Type",1);
150 con->setSolveParameters(list);
151 u_ptr->assign(nx, 1.0);
152 con->solve(*cp,*up,*zp,tol);
153 rnorm = cp->norm();
154 con->value(*cp,*up,*zp,tol);
155 cnorm = cp->norm();
156 errorFlag += ((cnorm > tol) ? 1 : 0) + ((rnorm > tol) ? 1 : 0);
157 *outStream << std::scientific << std::setprecision(8);
158 *outStream << std::endl;
159 *outStream << "Test SimOpt solve at feasible (u,z):" << std::endl;
160 *outStream << " Solver Residual = " << rnorm << std::endl;
161 *outStream << " ||c(u,z)|| = " << cnorm;
162 *outStream << std::endl << std::endl;
163 // Composite Step
164 list.sublist("SimOpt").sublist("Solve").set("Absolute Residual Tolerance",tol);
165 list.sublist("SimOpt").sublist("Solve").set("Solver Type",2);
166 con->setSolveParameters(list);
167 u_ptr->assign(nx, 1.0);
168 con->solve(*cp,*up,*zp,tol);
169 rnorm = cp->norm();
170 con->value(*cp,*up,*zp,tol);
171 cnorm = cp->norm();
172 tol *= 100.0; // Probably will not need this when we have composite step
173 errorFlag += ((cnorm > tol) ? 1 : 0) + ((rnorm > tol) ? 1 : 0);
174 *outStream << std::scientific << std::setprecision(8);
175 *outStream << std::endl;
176 *outStream << "Test SimOpt solve at feasible (u,z):" << std::endl;
177 *outStream << " Solver Residual = " << rnorm << std::endl;
178 *outStream << " ||c(u,z)|| = " << cnorm;
179 *outStream << std::endl << std::endl;
180 }
181 catch (std::logic_error& err) {
182 *outStream << err.what() << "\n";
183 errorFlag = -1000;
184 }; // end try
185
186 if (errorFlag != 0)
187 std::cout << "End Result: TEST FAILED\n";
188 else
189 std::cout << "End Result: TEST PASSED\n";
190
191 return 0;
192}
Defines a no-output stream class ROL::NullStream and a function makeStreamPtr which either wraps a re...
int main(int argc, char *argv[])
L2VectorPrimal< RealT > PrimalControlVector
H1VectorPrimal< RealT > DualConstraintVector
H1VectorPrimal< RealT > PrimalStateVector
H1VectorDual< RealT > DualStateVector
L2VectorDual< RealT > DualControlVector
H1VectorDual< RealT > PrimalConstraintVector
double RealT