47 #include "Epetra_MpiComm.h" 49 #include "Epetra_SerialComm.h" 53 #include "Ifpack2_Factory.hpp" 54 #include "BelosLinearProblem.hpp" 55 #include "BelosTpetraAdapter.hpp" 56 #include "BelosPseudoBlockCGSolMgr.hpp" 57 #include "BelosPseudoBlockGmresSolMgr.hpp" 58 #include "MatrixMarket_Tpetra.hpp" 64 #include "Teuchos_TimeMonitor.hpp" 68 typedef double MeshScalar;
69 typedef double BasisScalar;
72 typedef Tpetra::DefaultPlatform::DefaultPlatformType::NodeType
Node;
76 using Teuchos::ParameterList;
83 bool nonlinear_expansion =
false;
85 bool symmetric =
false;
89 MPI_Init(&argc,&
argv);
97 TEUCHOS_FUNC_TIME_MONITOR(
"Total PCE Calculation Time");
100 RCP<const Epetra_Comm> globalComm;
102 globalComm = rcp(
new Epetra_MpiComm(MPI_COMM_WORLD));
104 globalComm = rcp(
new Epetra_SerialComm);
106 MyPID = globalComm->MyPID();
109 Teuchos::Array< RCP<const Stokhos::OneDOrthogPolyBasis<LocalOrdinal,BasisScalar> > > bases(num_KL);
112 RCP<const Stokhos::CompletePolynomialBasis<LocalOrdinal,BasisScalar> > basis =
116 RCP<Stokhos::Sparse3Tensor<LocalOrdinal,BasisScalar> > Cijk;
117 if (nonlinear_expansion)
118 Cijk = basis->computeTripleProductTensor();
120 Cijk = basis->computeLinearTripleProductTensor();
121 RCP<Stokhos::OrthogPolyExpansion<LocalOrdinal,BasisScalar> > expansion =
125 std::cout <<
"Stochastic Galerkin expansion size = " << sz << std::endl;
130 num_spatial_procs = std::atoi(
argv[1]);
131 ParameterList parallelParams;
132 parallelParams.set(
"Number of Spatial Processors", num_spatial_procs);
137 RCP<Stokhos::ParallelData> sg_parallel_data =
140 RCP<const EpetraExt::MultiComm> sg_comm =
141 sg_parallel_data->getMultiComm();
142 RCP<const Epetra_Comm> app_comm =
143 sg_parallel_data->getSpatialComm();
146 RCP< Teuchos::Comm<int> > teuchos_app_comm;
148 RCP<const Epetra_MpiComm> app_mpi_comm =
149 Teuchos::rcp_dynamic_cast<
const Epetra_MpiComm>(app_comm);
150 RCP<const Teuchos::OpaqueWrapper<MPI_Comm> > raw_mpi_comm =
151 Teuchos::opaqueWrapper(app_mpi_comm->Comm());
152 teuchos_app_comm = rcp(
new Teuchos::MpiComm<int>(raw_mpi_comm));
154 teuchos_app_comm = rcp(
new Teuchos::SerialComm<int>());
159 RCP<problem_type> model =
160 rcp(
new problem_type(teuchos_app_comm, n, num_KL, s, mu,
161 nonlinear_expansion, symmetric));
165 typedef problem_type::Tpetra_Vector Tpetra_Vector;
166 typedef problem_type::Tpetra_CrsMatrix Tpetra_CrsMatrix;
167 typedef Tpetra::MatrixMarket::Writer<Tpetra_CrsMatrix> Writer;
168 RCP<const Tpetra_Vector> p = model->get_p_init(0);
169 RCP<Tpetra_Vector>
x = Tpetra::createVector<Scalar>(model->get_x_map());
171 RCP<Tpetra_Vector>
f = Tpetra::createVector<Scalar>(model->get_f_map());
172 RCP<Tpetra_Vector>
dx = Tpetra::createVector<Scalar>(model->get_x_map());
173 RCP<Tpetra_CrsMatrix> J = model->create_W();
176 Teuchos::ParameterList precParams;
177 std::string prec_name =
"RILUK";
178 precParams.set(
"fact: iluk level-of-fill", 1);
179 precParams.set(
"fact: iluk level-of-overlap", 0);
180 typedef Ifpack2::Preconditioner<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tprec;
181 Teuchos::RCP<Tprec> M;
182 Ifpack2::Factory factory;
183 M = factory.create<Tpetra_CrsMatrix>(prec_name, J);
184 M->setParameters(precParams);
187 model->computeResidual(*
x, *p, *
f);
188 model->computeJacobian(*
x, *p, *J);
194 f->norm2(Teuchos::arrayView(&norm_f,1));
196 std::cout <<
"\nInitial residual norm = " << norm_f << std::endl;
199 RCP<ParameterList> belosParams = rcp(
new ParameterList);
200 belosParams->set(
"Num Blocks", 20);
201 belosParams->set(
"Convergence Tolerance", 1e-12);
202 belosParams->set(
"Maximum Iterations", 1000);
203 belosParams->set(
"Verbosity", 33);
204 belosParams->set(
"Output Style", 1);
205 belosParams->set(
"Output Frequency", 1);
206 typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
208 typedef Belos::OperatorTraits<Scalar,MV,OP> BOPT;
209 typedef Belos::MultiVecTraits<Scalar,MV> BMVT;
210 typedef Belos::LinearProblem<Scalar,MV,OP> BLinProb;
211 RCP< BLinProb > problem = rcp(
new BLinProb(J,
dx,
f));
212 problem->setRightPrec(M);
213 problem->setProblem();
214 RCP<Belos::SolverManager<Scalar,MV,OP> > solver;
216 solver = rcp(
new Belos::PseudoBlockCGSolMgr<Scalar,MV,OP>(problem,
219 solver = rcp(
new Belos::PseudoBlockGmresSolMgr<Scalar,MV,OP>(problem,
223 Belos::ReturnType ret = solver->solve();
226 x->update(-1.0, *
dx, 1.0);
229 RCP<Tpetra_Vector>
g = Tpetra::createVector<Scalar>(model->get_g_map(0));
231 model->computeResidual(*
x, *p, *
f);
232 model->computeResponse(*
x, *p, *
g);
235 f->norm2(Teuchos::arrayView(&norm_f,1));
237 std::cout <<
"\nFinal residual norm = " << norm_f << std::endl;
240 std::cout <<
"\nResponse = " << std::endl;
241 Writer::writeDense(std::cout,
g);
245 Teuchos::TimeMonitor::summarize(std::cout);
246 Teuchos::TimeMonitor::zeroOutTimers();
250 catch (std::exception& e) {
251 std::cout << e.what() << std::endl;
254 std::cout << s << std::endl;
257 std::cout << s << std::endl;
260 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.