53 #include "Ifpack2_Factory.hpp" 54 #include "BelosLinearProblem.hpp" 55 #include "BelosTpetraAdapter.hpp" 56 #include "BelosPseudoBlockCGSolMgr.hpp" 57 #include "BelosPseudoBlockGmresSolMgr.hpp" 58 #include "Tpetra_Core.hpp" 59 #include "MatrixMarket_Tpetra.hpp" 65 #include "Teuchos_TimeMonitor.hpp" 69 typedef double MeshScalar;
70 typedef double BasisScalar;
77 using Teuchos::ParameterList;
84 bool nonlinear_expansion =
false;
86 bool symmetric =
false;
90 MPI_Init(&argc,&
argv);
98 TEUCHOS_FUNC_TIME_MONITOR(
"Total PCE Calculation Time");
101 RCP<const Epetra_Comm> globalComm;
107 MyPID = globalComm->MyPID();
110 Teuchos::Array< RCP<const Stokhos::OneDOrthogPolyBasis<LocalOrdinal,BasisScalar> > > bases(num_KL);
113 RCP<const Stokhos::CompletePolynomialBasis<LocalOrdinal,BasisScalar> > basis =
117 RCP<Stokhos::Sparse3Tensor<LocalOrdinal,BasisScalar> > Cijk;
118 if (nonlinear_expansion)
119 Cijk = basis->computeTripleProductTensor();
121 Cijk = basis->computeLinearTripleProductTensor();
122 RCP<Stokhos::OrthogPolyExpansion<LocalOrdinal,BasisScalar> > expansion =
126 std::cout <<
"Stochastic Galerkin expansion size = " << sz << std::endl;
131 num_spatial_procs = std::atoi(
argv[1]);
132 ParameterList parallelParams;
133 parallelParams.set(
"Number of Spatial Processors", num_spatial_procs);
138 RCP<Stokhos::ParallelData> sg_parallel_data =
141 RCP<const EpetraExt::MultiComm> sg_comm =
142 sg_parallel_data->getMultiComm();
143 RCP<const Epetra_Comm> app_comm =
144 sg_parallel_data->getSpatialComm();
147 RCP< Teuchos::Comm<int> > teuchos_app_comm;
149 RCP<const Epetra_MpiComm> app_mpi_comm =
151 RCP<const Teuchos::OpaqueWrapper<MPI_Comm> > raw_mpi_comm =
152 Teuchos::opaqueWrapper(app_mpi_comm->Comm());
153 teuchos_app_comm = rcp(
new Teuchos::MpiComm<int>(raw_mpi_comm));
155 teuchos_app_comm = rcp(
new Teuchos::SerialComm<int>());
160 RCP<problem_type> model =
161 rcp(
new problem_type(teuchos_app_comm, n, num_KL, s, mu,
162 nonlinear_expansion, symmetric));
166 typedef problem_type::Tpetra_Vector Tpetra_Vector;
167 typedef problem_type::Tpetra_CrsMatrix Tpetra_CrsMatrix;
168 typedef Tpetra::MatrixMarket::Writer<Tpetra_CrsMatrix> Writer;
169 RCP<const Tpetra_Vector> p = model->get_p_init(0);
170 RCP<Tpetra_Vector>
x = Tpetra::createVector<Scalar>(model->get_x_map());
172 RCP<Tpetra_Vector>
f = Tpetra::createVector<Scalar>(model->get_f_map());
173 RCP<Tpetra_Vector>
dx = Tpetra::createVector<Scalar>(model->get_x_map());
174 RCP<Tpetra_CrsMatrix> J = model->create_W();
177 Teuchos::ParameterList precParams;
178 std::string prec_name =
"RILUK";
179 precParams.set(
"fact: iluk level-of-fill", 1);
180 precParams.set(
"fact: iluk level-of-overlap", 0);
181 typedef Ifpack2::Preconditioner<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tprec;
182 Teuchos::RCP<Tprec> M;
183 Ifpack2::Factory factory;
184 M = factory.create<Tpetra_CrsMatrix>(prec_name, J);
185 M->setParameters(precParams);
188 model->computeResidual(*
x, *p, *
f);
189 model->computeJacobian(*
x, *p, *J);
195 f->norm2(Teuchos::arrayView(&norm_f,1));
197 std::cout <<
"\nInitial residual norm = " << norm_f << std::endl;
200 RCP<ParameterList> belosParams = rcp(
new ParameterList);
201 belosParams->set(
"Num Blocks", 20);
202 belosParams->set(
"Convergence Tolerance", 1e-12);
203 belosParams->set(
"Maximum Iterations", 1000);
204 belosParams->set(
"Verbosity", 33);
205 belosParams->set(
"Output Style", 1);
206 belosParams->set(
"Output Frequency", 1);
207 typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
208 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> MV;
209 typedef Belos::OperatorTraits<Scalar,MV,OP> BOPT;
210 typedef Belos::MultiVecTraits<Scalar,MV> BMVT;
211 typedef Belos::LinearProblem<Scalar,MV,OP> BLinProb;
212 RCP< BLinProb > problem = rcp(
new BLinProb(J,
dx,
f));
213 problem->setRightPrec(M);
214 problem->setProblem();
215 RCP<Belos::SolverManager<Scalar,MV,OP> > solver;
217 solver = rcp(
new Belos::PseudoBlockCGSolMgr<Scalar,MV,OP>(problem,
220 solver = rcp(
new Belos::PseudoBlockGmresSolMgr<Scalar,MV,OP>(problem,
224 Belos::ReturnType ret = solver->solve();
227 x->update(-1.0, *
dx, 1.0);
230 RCP<Tpetra_Vector>
g = Tpetra::createVector<Scalar>(model->get_g_map(0));
232 model->computeResidual(*
x, *p, *
f);
233 model->computeResponse(*
x, *p, *
g);
236 f->norm2(Teuchos::arrayView(&norm_f,1));
238 std::cout <<
"\nFinal residual norm = " << norm_f << std::endl;
241 std::cout <<
"\nResponse = " << std::endl;
242 Writer::writeDense(std::cout,
g);
246 Teuchos::TimeMonitor::summarize(std::cout);
247 Teuchos::TimeMonitor::zeroOutTimers();
251 catch (std::exception& e) {
252 std::cout << e.what() << std::endl;
255 std::cout << s << std::endl;
258 std::cout << s << std::endl;
261 std::cout <<
"Caught unknown exception!" <<std:: endl;
int main(int argc, char *argv[])
Orthogonal polynomial expansions limited to algebraic operations.
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType * x
KokkosClassic::DefaultNode::DefaultNodeType Node
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
Legendre polynomial basis.
ScalarType f(const Teuchos::Array< ScalarType > &x, double a, double b)
ScalarType g(const Teuchos::Array< ScalarType > &x, const ScalarType &y)
A linear 2-D diffusion problem.