Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
nox_example.cpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Stokhos Package
5// Copyright (2009) 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 Eric T. Phipps (etphipp@sandia.gov).
38//
39// ***********************************************************************
40// @HEADER
41
42#include <iostream>
43#include <sstream>
44
45// NOX
46#include "NOX.H"
47#include "NOX_Epetra.H"
48
49// Epetra communicator
50#ifdef HAVE_MPI
51#include "Epetra_MpiComm.h"
52#else
53#include "Epetra_SerialComm.h"
54#endif
55
56// Stokhos Stochastic Galerkin
57#include "Stokhos_Epetra.hpp"
58
59// Timing utilities
60#include "Teuchos_TimeMonitor.hpp"
61
62// Our model
63#include "SimpleME.hpp"
64
65int main(int argc, char *argv[]) {
66
67// Initialize MPI
68#ifdef HAVE_MPI
69 MPI_Init(&argc,&argv);
70#endif
71
72 int MyPID;
73
74 try {
75
76 // Create a communicator for Epetra objects
77 Teuchos::RCP<const Epetra_Comm> globalComm;
78#ifdef HAVE_MPI
79 globalComm = Teuchos::rcp(new Epetra_MpiComm(MPI_COMM_WORLD));
80#else
81 globalComm = Teuchos::rcp(new Epetra_SerialComm);
82#endif
83 MyPID = globalComm->MyPID();
84
85 // Create Stochastic Galerkin basis and expansion
86 int p = 5;
87 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > bases(1);
88 bases[0] = Teuchos::rcp(new Stokhos::LegendreBasis<int,double>(p));
89 Teuchos::RCP<const Stokhos::CompletePolynomialBasis<int,double> > basis =
90 Teuchos::rcp(new Stokhos::CompletePolynomialBasis<int,double>(bases));
91 int sz = basis->size();
92 Teuchos::RCP<Stokhos::Sparse3Tensor<int,double> > Cijk;
93 Cijk = basis->computeTripleProductTensor();
94 Teuchos::RCP<Stokhos::OrthogPolyExpansion<int,double> > expansion =
96 Cijk));
97 if (MyPID == 0)
98 std::cout << "Stochastic Galerkin expansion size = " << sz << std::endl;
99
100 // Create stochastic parallel distribution
101 int num_spatial_procs = -1;
102 Teuchos::ParameterList parallelParams;
103 parallelParams.set("Number of Spatial Processors", num_spatial_procs);
104 Teuchos::RCP<Stokhos::ParallelData> sg_parallel_data =
105 Teuchos::rcp(new Stokhos::ParallelData(basis, Cijk, globalComm,
106 parallelParams));
107 Teuchos::RCP<const EpetraExt::MultiComm> sg_comm =
108 sg_parallel_data->getMultiComm();
109 Teuchos::RCP<const Epetra_Comm> app_comm =
110 sg_parallel_data->getSpatialComm();
111
112 // Create application model evaluator
113 Teuchos::RCP<EpetraExt::ModelEvaluator> model =
114 Teuchos::rcp(new SimpleME(app_comm));
115
116 // Setup stochastic Galerkin algorithmic parameters
117 Teuchos::RCP<Teuchos::ParameterList> sgParams =
118 Teuchos::rcp(new Teuchos::ParameterList);
119 sgParams->set("Jacobian Method", "Matrix Free");
120 sgParams->set("Mean Preconditioner Type", "Ifpack");
121
122 // Create stochastic Galerkin model evaluator
123 Teuchos::RCP<Stokhos::SGModelEvaluator> sg_model =
124 Teuchos::rcp(new Stokhos::SGModelEvaluator(model, basis, Teuchos::null,
125 expansion, sg_parallel_data,
126 sgParams));
127
128 // Stochastic Galerkin initial guess
129 // Set the mean to the deterministic initial guess, higher-order terms
130 // to zero
131 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> x_init_sg =
132 sg_model->create_x_sg();
133 x_init_sg->init(0.0);
134 (*x_init_sg)[0] = *(model->get_x_init());
135 sg_model->set_x_sg_init(*x_init_sg);
136
137 // Stochastic Galerkin parameters
138 // Linear expansion with the mean given by the deterministic initial
139 // parameter values, linear terms equal to 1, and higher order terms
140 // equal to zero.
141 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> p_init_sg =
142 sg_model->create_p_sg(0);
143 p_init_sg->init(0.0);
144 (*p_init_sg)[0] = *(model->get_p_init(0));
145 for (int i=0; i<model->get_p_map(0)->NumMyElements(); i++)
146 (*p_init_sg)[i+1][i] = 1.0;
147 sg_model->set_p_sg_init(0, *p_init_sg);
148 std::cout << "Stochatic Galerkin parameter expansion = " << std::endl
149 << *p_init_sg << std::endl;
150
151 // Set up NOX parameters
152 Teuchos::RCP<Teuchos::ParameterList> noxParams =
153 Teuchos::rcp(new Teuchos::ParameterList);
154
155 // Set the nonlinear solver method
156 noxParams->set("Nonlinear Solver", "Line Search Based");
157
158 // Set the printing parameters in the "Printing" sublist
159 Teuchos::ParameterList& printParams = noxParams->sublist("Printing");
160 printParams.set("MyPID", MyPID);
161 printParams.set("Output Precision", 3);
162 printParams.set("Output Processor", 0);
163 printParams.set("Output Information",
164 NOX::Utils::OuterIteration +
165 NOX::Utils::OuterIterationStatusTest +
166 NOX::Utils::InnerIteration +
167 //NOX::Utils::Parameters +
168 //NOX::Utils::Details +
169 NOX::Utils::LinearSolverDetails +
170 NOX::Utils::Warning +
171 NOX::Utils::Error);
172
173 // Create printing utilities
174 NOX::Utils utils(printParams);
175
176 // Sublist for line search
177 Teuchos::ParameterList& searchParams = noxParams->sublist("Line Search");
178 searchParams.set("Method", "Full Step");
179
180 // Sublist for direction
181 Teuchos::ParameterList& dirParams = noxParams->sublist("Direction");
182 dirParams.set("Method", "Newton");
183 Teuchos::ParameterList& newtonParams = dirParams.sublist("Newton");
184 newtonParams.set("Forcing Term Method", "Constant");
185
186 // Sublist for linear solver for the Newton method
187 Teuchos::ParameterList& lsParams = newtonParams.sublist("Linear Solver");
188 lsParams.set("Aztec Solver", "GMRES");
189 lsParams.set("Max Iterations", 100);
190 lsParams.set("Size of Krylov Subspace", 100);
191 lsParams.set("Tolerance", 1e-4);
192 lsParams.set("Output Frequency", 10);
193 lsParams.set("Preconditioner", "Ifpack");
194
195 // Sublist for convergence tests
196 Teuchos::ParameterList& statusParams = noxParams->sublist("Status Tests");
197 statusParams.set("Test Type", "Combo");
198 statusParams.set("Number of Tests", 2);
199 statusParams.set("Combo Type", "OR");
200 Teuchos::ParameterList& normF = statusParams.sublist("Test 0");
201 normF.set("Test Type", "NormF");
202 normF.set("Tolerance", 1e-10);
203 normF.set("Scale Type", "Scaled");
204 Teuchos::ParameterList& maxIters = statusParams.sublist("Test 1");
205 maxIters.set("Test Type", "MaxIters");
206 maxIters.set("Maximum Iterations", 10);
207
208 // Create NOX interface
209 Teuchos::RCP<NOX::Epetra::ModelEvaluatorInterface> nox_interface =
210 Teuchos::rcp(new NOX::Epetra::ModelEvaluatorInterface(sg_model));
211
212 // Create NOX linear system object
213 Teuchos::RCP<const Epetra_Vector> u = sg_model->get_x_init();
214 Teuchos::RCP<Epetra_Operator> A = sg_model->create_W();
215 Teuchos::RCP<NOX::Epetra::Interface::Required> iReq = nox_interface;
216 Teuchos::RCP<NOX::Epetra::Interface::Jacobian> iJac = nox_interface;
217 Teuchos::RCP<NOX::Epetra::LinearSystemAztecOO> linsys;
218 Teuchos::RCP<Epetra_Operator> M = sg_model->create_WPrec()->PrecOp;
219 Teuchos::RCP<NOX::Epetra::Interface::Preconditioner> iPrec = nox_interface;
220 lsParams.set("Preconditioner", "User Defined");
221 linsys =
222 Teuchos::rcp(new NOX::Epetra::LinearSystemAztecOO(printParams, lsParams,
223 iJac, A, iPrec, M,
224 *u));
225 // linsys =
226 // Teuchos::rcp(new NOX::Epetra::LinearSystemAztecOO(printParams, lsParams,
227 // iReq, iJac, A, *u));
228
229 // Build NOX group
230 Teuchos::RCP<NOX::Epetra::Group> grp =
231 Teuchos::rcp(new NOX::Epetra::Group(printParams, iReq, *u, linsys));
232
233 // Create the Solver convergence test
234 Teuchos::RCP<NOX::StatusTest::Generic> statusTests =
235 NOX::StatusTest::buildStatusTests(statusParams, utils);
236
237 // Create the solver
238 Teuchos::RCP<NOX::Solver::Generic> solver =
239 NOX::Solver::buildSolver(grp, statusTests, noxParams);
240
241 // Solve the system
242 NOX::StatusTest::StatusType status = solver->solve();
243
244 // Get final solution
245 const NOX::Epetra::Group& finalGroup =
246 dynamic_cast<const NOX::Epetra::Group&>(solver->getSolutionGroup());
247 const Epetra_Vector& finalSolution =
248 (dynamic_cast<const NOX::Epetra::Vector&>(finalGroup.getX())).getEpetraVector();
249
250 // Convert block Epetra_Vector to orthogonal polynomial of Epetra_Vector's
251 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> x_sg =
252 sg_model->create_x_sg(View, &finalSolution);
253
254 utils.out() << "Final Solution (block vector) = " << std::endl;
255 std::cout << finalSolution << std::endl;
256 utils.out() << "Final Solution (polynomial) = " << std::endl;
257 std::cout << *x_sg << std::endl;
258
259 if (status == NOX::StatusTest::Converged && MyPID == 0)
260 utils.out() << "Example Passed!" << std::endl;
261
262 }
263
264 catch (std::exception& e) {
265 std::cout << e.what() << std::endl;
266 }
267 catch (std::string& s) {
268 std::cout << s << std::endl;
269 }
270 catch (char *s) {
271 std::cout << s << std::endl;
272 }
273 catch (...) {
274 std::cout << "Caught unknown exception!" << std::endl;
275 }
276
277#ifdef HAVE_MPI
278 MPI_Finalize() ;
279#endif
280
281}
Orthogonal polynomial expansions limited to algebraic operations.
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
Legendre polynomial basis.
Nonlinear, stochastic Galerkin ModelEvaluator.
int main(int argc, char *argv[])