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