47 #include "NOX_Epetra.H" 51 #include "Epetra_MpiComm.h" 53 #include "Epetra_SerialComm.h" 60 #include "Teuchos_TimeMonitor.hpp" 69 MPI_Init(&argc,&
argv);
77 Teuchos::RCP<const Epetra_Comm> globalComm;
79 globalComm = Teuchos::rcp(
new Epetra_MpiComm(MPI_COMM_WORLD));
81 globalComm = Teuchos::rcp(
new Epetra_SerialComm);
83 MyPID = globalComm->MyPID();
87 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > bases(1);
89 Teuchos::RCP<const Stokhos::CompletePolynomialBasis<int,double> > basis =
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 =
98 std::cout <<
"Stochastic Galerkin expansion size = " << sz << std::endl;
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 =
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();
113 Teuchos::RCP<EpetraExt::ModelEvaluator> model =
114 Teuchos::rcp(
new SimpleME(app_comm));
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");
123 Teuchos::RCP<Stokhos::SGModelEvaluator> sg_model =
125 expansion, sg_parallel_data,
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);
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;
152 Teuchos::RCP<Teuchos::ParameterList> noxParams =
153 Teuchos::rcp(
new Teuchos::ParameterList);
156 noxParams->set(
"Nonlinear Solver",
"Line Search Based");
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 +
169 NOX::Utils::LinearSolverDetails +
170 NOX::Utils::Warning +
174 NOX::Utils utils(printParams);
177 Teuchos::ParameterList& searchParams = noxParams->sublist(
"Line Search");
178 searchParams.set(
"Method",
"Full Step");
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");
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");
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);
209 Teuchos::RCP<NOX::Epetra::ModelEvaluatorInterface> nox_interface =
210 Teuchos::rcp(
new NOX::Epetra::ModelEvaluatorInterface(sg_model));
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");
222 Teuchos::rcp(
new NOX::Epetra::LinearSystemAztecOO(printParams, lsParams,
230 Teuchos::RCP<NOX::Epetra::Group> grp =
231 Teuchos::rcp(
new NOX::Epetra::Group(printParams, iReq, *u, linsys));
234 Teuchos::RCP<NOX::StatusTest::Generic> statusTests =
235 NOX::StatusTest::buildStatusTests(statusParams, utils);
238 Teuchos::RCP<NOX::Solver::Generic> solver =
239 NOX::Solver::buildSolver(grp, statusTests, noxParams);
242 NOX::StatusTest::StatusType status = solver->solve();
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();
251 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> x_sg =
252 sg_model->create_x_sg(View, &finalSolution);
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;
259 if (status == NOX::StatusTest::Converged && MyPID == 0)
260 utils.out() <<
"Example Passed!" << std::endl;
264 catch (std::exception& e) {
265 std::cout << e.what() << std::endl;
267 catch (std::string& s) {
268 std::cout << s << std::endl;
271 std::cout << s << std::endl;
274 std::cout <<
"Caught unknown exception!" << std::endl;
Orthogonal polynomial expansions limited to algebraic operations.
Nonlinear, stochastic Galerkin ModelEvaluator.
int main(int argc, char *argv[])
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
Legendre polynomial basis.