ROL
function/test_04.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 
14 #include "ROL_Stream.hpp"
15 
16 #include "Teuchos_GlobalMPISession.hpp"
17 #include "Teuchos_Comm.hpp"
18 #include "Teuchos_DefaultComm.hpp"
19 #include "Teuchos_CommHelpers.hpp"
20 
21 #include <iostream>
22 #include <fstream>
23 #include <algorithm>
24 
25 #include "test_04.hpp"
26 
27 typedef double RealT;
34 
35 int main(int argc, char *argv[]) {
36 
37  Teuchos::GlobalMPISession mpiSession(&argc, &argv);
38 
39  // This little trick lets us print to std::cout only if a (dummy) command-line argument is provided.
40  int iprint = argc - 1;
41  bool print = (iprint>0);
42  ROL::Ptr<std::ostream> outStream;
43  ROL::nullstream bhs; // outputs nothing
44  if (print)
45  outStream = ROL::makePtrFromRef(std::cout);
46  else
47  outStream = ROL::makePtrFromRef(bhs);
48 
49  int errorFlag = 0;
50 
51  // *** Example body.
52 
53  try {
54  /*************************************************************************/
55  /************* INITIALIZE BURGERS FEM CLASS ******************************/
56  /*************************************************************************/
57  int nx = 512; // Set spatial discretization.
58  RealT nl = 1.0; // Nonlinearity parameter (1 = Burgers, 0 = linear).
59  RealT cH1 = 1.0; // Scale for derivative term in H1 norm.
60  RealT cL2 = 0.0; // Scale for mass term in H1 norm.
61  ROL::Ptr<BurgersFEM<RealT> > fem
62  = ROL::makePtr<BurgersFEM<RealT>>(nx,nl,cH1,cL2);
63  fem->test_inverse_mass(*outStream);
64  fem->test_inverse_H1(*outStream);
65  /*************************************************************************/
66  /************* INITIALIZE SIMOPT CONSTRAINT ******************************/
67  /*************************************************************************/
68  bool hess = true;
69  ROL::Ptr<ROL::Constraint_SimOpt<RealT> > con
70  = ROL::makePtr<Constraint_BurgersControl<RealT>>(fem,hess);
71  /*************************************************************************/
72  /************* INITIALIZE VECTOR STORAGE *********************************/
73  /*************************************************************************/
74  // INITIALIZE CONTROL VECTORS
75  ROL::Ptr<std::vector<RealT> > z_ptr
76  = ROL::makePtr<std::vector<RealT>>(nx+2, 0.0);
77  ROL::Ptr<ROL::Vector<RealT> > zp
78  = ROL::makePtr<PrimalControlVector>(z_ptr,fem);
79  // INITIALIZE STATE VECTORS
80  ROL::Ptr<std::vector<RealT> > u_ptr
81  = ROL::makePtr<std::vector<RealT>>(nx, 1.0);
82  ROL::Ptr<ROL::Vector<RealT> > up
83  = ROL::makePtr<PrimalStateVector>(u_ptr,fem);
84  // INITIALIZE CONSTRAINT VECTORS
85  ROL::Ptr<std::vector<RealT> > c_ptr
86  = ROL::makePtr<std::vector<RealT>>(nx, 1.0);
87  ROL::Ptr<ROL::Vector<RealT> > cp
88  = ROL::makePtr<PrimalConstraintVector>(c_ptr,fem);
89  /*************************************************************************/
90  /************* CHECK DERIVATIVES AND CONSISTENCY *************************/
91  /*************************************************************************/
92  RealT rnorm(0), cnorm(0);
93  ROL::ParameterList list;
94  RealT tol = std::sqrt(ROL::ROL_EPSILON<RealT>());
95  list.sublist("SimOpt").sublist("Solve").set("Output Iteration History",print);
96  list.sublist("SimOpt").sublist("Solve").set("Step Tolerance",ROL::ROL_EPSILON<RealT>());
97  // Newton
98  list.sublist("SimOpt").sublist("Solve").set("Absolute Residual Tolerance",tol);
99  list.sublist("SimOpt").sublist("Solve").set("Solver Type",0);
100  con->setSolveParameters(list);
101  u_ptr->assign(nx, 1.0);
102  con->solve(*cp,*up,*zp,tol);
103  rnorm = cp->norm();
104  con->value(*cp,*up,*zp,tol);
105  cnorm = cp->norm();
106  errorFlag += ((cnorm > tol) ? 1 : 0) + ((rnorm > tol) ? 1 : 0);
107  *outStream << std::scientific << std::setprecision(8);
108  *outStream << std::endl;
109  *outStream << "Test SimOpt solve at feasible (u,z):" << std::endl;
110  *outStream << " Solver Residual = " << rnorm << std::endl;
111  *outStream << " ||c(u,z)|| = " << cnorm;
112  *outStream << std::endl << std::endl;
113  // Levenberg-Marquardt
114  list.sublist("SimOpt").sublist("Solve").set("Absolute Residual Tolerance",1e-4*tol);
115  list.sublist("SimOpt").sublist("Solve").set("Solver Type",1);
116  con->setSolveParameters(list);
117  u_ptr->assign(nx, 1.0);
118  con->solve(*cp,*up,*zp,tol);
119  rnorm = cp->norm();
120  con->value(*cp,*up,*zp,tol);
121  cnorm = cp->norm();
122  errorFlag += ((cnorm > tol) ? 1 : 0) + ((rnorm > tol) ? 1 : 0);
123  *outStream << std::scientific << std::setprecision(8);
124  *outStream << std::endl;
125  *outStream << "Test SimOpt solve at feasible (u,z):" << std::endl;
126  *outStream << " Solver Residual = " << rnorm << std::endl;
127  *outStream << " ||c(u,z)|| = " << cnorm;
128  *outStream << std::endl << std::endl;
129  // Composite Step
130  list.sublist("SimOpt").sublist("Solve").set("Absolute Residual Tolerance",tol);
131  list.sublist("SimOpt").sublist("Solve").set("Solver Type",2);
132  con->setSolveParameters(list);
133  u_ptr->assign(nx, 1.0);
134  con->solve(*cp,*up,*zp,tol);
135  rnorm = cp->norm();
136  con->value(*cp,*up,*zp,tol);
137  cnorm = cp->norm();
138  tol *= 100.0; // Probably will not need this when we have composite step
139  errorFlag += ((cnorm > tol) ? 1 : 0) + ((rnorm > tol) ? 1 : 0);
140  *outStream << std::scientific << std::setprecision(8);
141  *outStream << std::endl;
142  *outStream << "Test SimOpt solve at feasible (u,z):" << std::endl;
143  *outStream << " Solver Residual = " << rnorm << std::endl;
144  *outStream << " ||c(u,z)|| = " << cnorm;
145  *outStream << std::endl << std::endl;
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 }
L2VectorPrimal< RealT > PrimalControlVector
Defines a no-output stream class ROL::NullStream and a function makeStreamPtr which either wraps a re...
H1VectorDual< RealT > DualStateVector
L2VectorDual< RealT > DualControlVector
H1VectorDual< RealT > PrimalConstraintVector
basic_nullstream< char, char_traits< char >> nullstream
Definition: ROL_Stream.hpp:38
int main(int argc, char *argv[])
H1VectorPrimal< RealT > PrimalStateVector
H1VectorPrimal< RealT > DualConstraintVector