diff --git a/ippl/test/AMR/CMakeLists.txt b/ippl/test/AMR/CMakeLists.txt
index a88a076abda5f7acb0f337148d97f8cbb7bf4379..47655f299e84cfa79678e499f77624a2ebc22561 100644
--- a/ippl/test/AMR/CMakeLists.txt
+++ b/ippl/test/AMR/CMakeLists.txt
@@ -50,9 +50,9 @@ IF ( ENABLE_AMR )
add_library(multigrid
trilinos/AmrMultiGridCore.h
trilinos/AmrMultiGridDefs.h
- trilinos/BelosBottomSolver.cpp
- trilinos/AmesosBottomSolver.cpp
- trilinos/MueLuBottomSolver.cpp
+ trilinos/BelosBottomSolver.hpp
+ trilinos/AmesosBottomSolver.hpp
+ trilinos/MueLuBottomSolver.hpp
# trilinos/FFTBottomSolver.cpp
trilinos/AmrSmoother.cpp
trilinos/AmrMultiGrid.cpp
diff --git a/ippl/test/AMR/trilinos/AmesosBottomSolver.cpp b/ippl/test/AMR/trilinos/AmesosBottomSolver.cpp
deleted file mode 100644
index 94c7705620a86897845047f1f4d6efae8734d311..0000000000000000000000000000000000000000
--- a/ippl/test/AMR/trilinos/AmesosBottomSolver.cpp
+++ /dev/null
@@ -1,38 +0,0 @@
-#include "AmesosBottomSolver.h"
-
-AmesosBottomSolver::AmesosBottomSolver(std::string solvertype)
- : solvertype_m(solvertype)
-{ }
-
-
-AmesosBottomSolver::~AmesosBottomSolver() {
- solver_mp = Teuchos::null;
-}
-
-
-void AmesosBottomSolver::solve(const Teuchos::RCP<mv_t>& x,
- const Teuchos::RCP<mv_t>& b)
-{
- /*
- * solve linear system Ax = b
- */
- solver_mp->solve(x.get(), b.get());
-}
-
-
-void AmesosBottomSolver::setOperator(const Teuchos::RCP<matrix_t>& A) {
- try {
- solver_mp = Amesos2::create<matrix_t, mv_t>(solvertype_m, A);
- } catch(const std::invalid_argument& ex) {
- //TODO change to gmsg IPPL when built-in in OPAL
- std::cerr << ex.what() << std::endl;
- }
-
- solver_mp->symbolicFactorization();
- solver_mp->numericFactorization();
-}
-
-
-std::size_t AmesosBottomSolver::getNumIters() {
- return 1; // direct solvers do only one step
-}
diff --git a/ippl/test/AMR/trilinos/AmesosBottomSolver.h b/ippl/test/AMR/trilinos/AmesosBottomSolver.h
index 636ae7a76c64a048ec93bd3aa17e7a956508c0b3..7872a16ddf1c3151a747fef1486d9345f8e45932 100644
--- a/ippl/test/AMR/trilinos/AmesosBottomSolver.h
+++ b/ippl/test/AMR/trilinos/AmesosBottomSolver.h
@@ -8,8 +8,10 @@
#include <string>
/// Interface to Amesos2 solvers of the Trilinos package
+template <class Level>
class AmesosBottomSolver : public BottomSolver<Teuchos::RCP<amr::matrix_t>,
- Teuchos::RCP<amr::multivector_t> >
+ Teuchos::RCP<amr::multivector_t>,
+ Level>
{
public:
typedef amr::matrix_t matrix_t;
@@ -30,7 +32,8 @@ public:
void solve(const Teuchos::RCP<mv_t>& x,
const Teuchos::RCP<mv_t>& b);
- void setOperator(const Teuchos::RCP<matrix_t>& A);
+ void setOperator(const Teuchos::RCP<matrix_t>& A,
+ Level* level_p = nullptr);
std::size_t getNumIters();
@@ -41,4 +44,6 @@ private:
Teuchos::RCP<solver_t> solver_mp; ///< solver instance
};
+#include "AmesosBottomSolver.hpp"
+
#endif
diff --git a/ippl/test/AMR/trilinos/AmesosBottomSolver.hpp b/ippl/test/AMR/trilinos/AmesosBottomSolver.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..743fa21f570dbc45f09f28336814e4488767392b
--- /dev/null
+++ b/ippl/test/AMR/trilinos/AmesosBottomSolver.hpp
@@ -0,0 +1,43 @@
+template <class Level>
+AmesosBottomSolver<Level>::AmesosBottomSolver(std::string solvertype)
+ : solvertype_m(solvertype)
+{ }
+
+
+template <class Level>
+AmesosBottomSolver<Level>::~AmesosBottomSolver() {
+ solver_mp = Teuchos::null;
+}
+
+
+template <class Level>
+void AmesosBottomSolver<Level>::solve(const Teuchos::RCP<mv_t>& x,
+ const Teuchos::RCP<mv_t>& b)
+{
+ /*
+ * solve linear system Ax = b
+ */
+ solver_mp->solve(x.get(), b.get());
+}
+
+
+template <class Level>
+void AmesosBottomSolver<Level>::setOperator(const Teuchos::RCP<matrix_t>& A,
+ Level* level_p)
+{
+ try {
+ solver_mp = Amesos2::create<matrix_t, mv_t>(solvertype_m, A);
+ } catch(const std::invalid_argument& ex) {
+ //TODO change to gmsg IPPL when built-in in OPAL
+ std::cerr << ex.what() << std::endl;
+ }
+
+ solver_mp->symbolicFactorization();
+ solver_mp->numericFactorization();
+}
+
+
+template <class Level>
+std::size_t AmesosBottomSolver<Level>::getNumIters() {
+ return 1; // direct solvers do only one step
+}
diff --git a/ippl/test/AMR/trilinos/AmrMultiGrid.cpp b/ippl/test/AMR/trilinos/AmrMultiGrid.cpp
index 8d8fe5b43a48b8b55ec8d976ff37f052c72261cd..f73606d7d6647fdd61f7fd08083dd2cb7e42aa0a 100644
--- a/ippl/test/AMR/trilinos/AmrMultiGrid.cpp
+++ b/ippl/test/AMR/trilinos/AmrMultiGrid.cpp
@@ -33,7 +33,7 @@ AmrMultiGrid::AmrMultiGrid(AmrOpal* itsAmrObject_p,
bIter_m(0),
maxiter_m(100),
nSweeps_m(nSweeps),
- mglevel_m(0),
+ mglevel_m(1),
lbase_m(0),
lfine_m(0),
nlevel_m(1),
@@ -49,6 +49,8 @@ AmrMultiGrid::AmrMultiGrid(AmrOpal* itsAmrObject_p,
this->initTimer_m();
#endif
+ mglevel_m[0].reset(nullptr);
+
const Boundary bcs[AMREX_SPACEDIM] = {
this->convertToEnumBoundary_m(bcx),
this->convertToEnumBoundary_m(bcy),
@@ -672,7 +674,7 @@ void AmrMultiGrid::setup_m(const amrex::Array<AmrField_u>& rho,
if ( isNewOperator ) {
// set the bottom solve operator
- solver_mp->setOperator(mglevel_m[lbase_m]->Anf_p);
+ solver_mp->setOperator(mglevel_m[lbase_m]->Anf_p, mglevel_m[0].get());
}
@@ -1853,62 +1855,62 @@ void AmrMultiGrid::initBaseSolver_m(const BaseSolver& solver)
switch ( solver ) {
// Belos solvers
case BaseSolver::BICGSTAB:
- solver_mp.reset( new BelosBottomSolver("BICGSTAB", prec_mp) );
+ solver_mp.reset( new BelosBottomSolver<AmrMultiGridLevel_t>("BICGSTAB", prec_mp) );
break;
case BaseSolver::MINRES:
- solver_mp.reset( new BelosBottomSolver("MINRES", prec_mp) );
+ solver_mp.reset( new BelosBottomSolver<AmrMultiGridLevel_t>("MINRES", prec_mp) );
break;
case BaseSolver::PCPG:
- solver_mp.reset( new BelosBottomSolver("PCPG", prec_mp) );
+ solver_mp.reset( new BelosBottomSolver<AmrMultiGridLevel_t>("PCPG", prec_mp) );
break;
case BaseSolver::CG:
- solver_mp.reset( new BelosBottomSolver("Pseudoblock CG", prec_mp) );
+ solver_mp.reset( new BelosBottomSolver<AmrMultiGridLevel_t>("Pseudoblock CG", prec_mp) );
break;
case BaseSolver::GMRES:
- solver_mp.reset( new BelosBottomSolver("Pseudoblock GMRES", prec_mp) );
+ solver_mp.reset( new BelosBottomSolver<AmrMultiGridLevel_t>("Pseudoblock GMRES", prec_mp) );
break;
case BaseSolver::STOCHASTIC_CG:
- solver_mp.reset( new BelosBottomSolver("Stochastic CG", prec_mp) );
+ solver_mp.reset( new BelosBottomSolver<AmrMultiGridLevel_t>("Stochastic CG", prec_mp) );
break;
case BaseSolver::RECYCLING_CG:
- solver_mp.reset( new BelosBottomSolver("RCG", prec_mp) );
+ solver_mp.reset( new BelosBottomSolver<AmrMultiGridLevel_t>("RCG", prec_mp) );
break;
case BaseSolver::RECYCLING_GMRES:
- solver_mp.reset( new BelosBottomSolver("GCRODR", prec_mp) );
+ solver_mp.reset( new BelosBottomSolver<AmrMultiGridLevel_t>("GCRODR", prec_mp) );
break;
// Amesos2 solvers
#ifdef HAVE_AMESOS2_KLU2
case BaseSolver::KLU2:
- solver_mp.reset( new AmesosBottomSolver("klu2") );
+ solver_mp.reset( new AmesosBottomSolver<AmrMultiGridLevel_t>("klu2") );
break;
#endif
#ifdef HAVE_AMESOS2_SUPERLU
case BaseSolver::SUPERLU:
- solver_mp.reset( new AmesosBottomSolver("superlu") );
+ solver_mp.reset( new AmesosBottomSolver<AmrMultiGridLevel_t>("superlu") );
break;
#endif
#ifdef HAVE_AMESOS2_UMFPACK
case BaseSolver::UMFPACK:
- solver_mp.reset( new AmesosBottomSolver("umfpack") );
+ solver_mp.reset( new AmesosBottomSolver<AmrMultiGridLevel_t>("umfpack") );
break;
#endif
#ifdef HAVE_AMESOS2_PARDISO_MKL
case BaseSolver::PARDISO_MKL:
- solver_mp.reset( new AmesosBottomSolver("pardiso_mkl") );
+ solver_mp.reset( new AmesosBottomSolver<AmrMultiGridLevel_t>("pardiso_mkl") );
break;
#endif
#ifdef HAVE_AMESOS2_MUMPS
case BaseSolver::MUMPS:
- solver_mp.reset( new AmesosBottomSolver("mumps") );
+ solver_mp.reset( new AmesosBottomSolver<AmrMultiGridLevel_t>("mumps") );
break;
#endif
#ifdef HAVE_AMESOS2_LAPACK
case BaseSolver::LAPACK:
- solver_mp.reset( new AmesosBottomSolver("lapack") );
+ solver_mp.reset( new AmesosBottomSolver<AmrMultiGridLevel_t>("lapack") );
break;
#endif
case BaseSolver::SA:
- solver_mp.reset( new MueLuBottomSolver() );
+ solver_mp.reset( new MueLuBottomSolver<AmrMultiGridLevel_t>() );
break;
default:
throw OpalException("AmrMultiGrid::initBaseSolver_m()",
diff --git a/ippl/test/AMR/trilinos/AmrMultiGrid.h b/ippl/test/AMR/trilinos/AmrMultiGrid.h
index c76f67287b9a73c8efb5ceefd9db87bc87f7b2db..9113b77d3e7016445c1b5639befc1ccfa486057a 100644
--- a/ippl/test/AMR/trilinos/AmrMultiGrid.h
+++ b/ippl/test/AMR/trilinos/AmrMultiGrid.h
@@ -48,7 +48,8 @@ public:
typedef BottomSolver<
Teuchos::RCP<matrix_t>,
- Teuchos::RCP<mv_t>
+ Teuchos::RCP<mv_t>,
+ AmrMultiGridLevel_t
> bsolver_t;
typedef AmrPreconditioner<matrix_t> preconditioner_t;
diff --git a/ippl/test/AMR/trilinos/BelosBottomSolver.h b/ippl/test/AMR/trilinos/BelosBottomSolver.h
index 17539f9d0b6180440d9c08ef3f7917dea83de9eb..c12537ad0d6b4849f5288214528d72c845b6cb63 100644
--- a/ippl/test/AMR/trilinos/BelosBottomSolver.h
+++ b/ippl/test/AMR/trilinos/BelosBottomSolver.h
@@ -13,8 +13,10 @@
#include <string>
/// Interface to Belos solvers of the Trilinos package
+template <class Level>
class BelosBottomSolver : public BottomSolver<Teuchos::RCP<amr::matrix_t>,
- Teuchos::RCP<amr::multivector_t> >
+ Teuchos::RCP<amr::multivector_t>,
+ Level>
{
public:
typedef amr::matrix_t matrix_t;
@@ -45,7 +47,8 @@ public:
void solve(const Teuchos::RCP<mv_t>& x,
const Teuchos::RCP<mv_t>& b);
- void setOperator(const Teuchos::RCP<matrix_t>& A);
+ void setOperator(const Teuchos::RCP<matrix_t>& A,
+ Level* level_p = nullptr);
std::size_t getNumIters();
@@ -68,4 +71,6 @@ private:
int maxiter_m;
};
+#include "BelosBottomSolver.hpp"
+
#endif
diff --git a/ippl/test/AMR/trilinos/BelosBottomSolver.cpp b/ippl/test/AMR/trilinos/BelosBottomSolver.hpp
similarity index 76%
rename from ippl/test/AMR/trilinos/BelosBottomSolver.cpp
rename to ippl/test/AMR/trilinos/BelosBottomSolver.hpp
index 3799a354045548a8260c26868b18a12dffba4e3a..1b465e5f7c0da8bba5204b8ed5b3772a568b3c58 100644
--- a/ippl/test/AMR/trilinos/BelosBottomSolver.cpp
+++ b/ippl/test/AMR/trilinos/BelosBottomSolver.hpp
@@ -1,10 +1,9 @@
-#include "BelosBottomSolver.h"
-
#include "Ippl.h"
#include "Utilities/OpalException.h"
-BelosBottomSolver::BelosBottomSolver(std::string solvertype,
- const std::shared_ptr<prec_t>& prec_p)
+template <class Level>
+BelosBottomSolver<Level>::BelosBottomSolver(std::string solvertype,
+ const std::shared_ptr<prec_t>& prec_p)
: problem_mp( Teuchos::rcp( new problem_t() ) ),
prec_mp(prec_p),
reltol_m(1.0e-9),
@@ -14,15 +13,17 @@ BelosBottomSolver::BelosBottomSolver(std::string solvertype,
}
-BelosBottomSolver::~BelosBottomSolver() {
+template <class Level>
+BelosBottomSolver<Level>::~BelosBottomSolver() {
problem_mp = Teuchos::null;
params_mp = Teuchos::null;
solver_mp = Teuchos::null;
}
-void BelosBottomSolver::solve(const Teuchos::RCP<mv_t>& x,
- const Teuchos::RCP<mv_t>& b)
+template <class Level>
+void BelosBottomSolver<Level>::solve(const Teuchos::RCP<mv_t>& x,
+ const Teuchos::RCP<mv_t>& b)
{
/*
* solve linear system Ax = b
@@ -49,7 +50,10 @@ void BelosBottomSolver::solve(const Teuchos::RCP<mv_t>& x,
}
-void BelosBottomSolver::setOperator(const Teuchos::RCP<matrix_t>& A) {
+template <class Level>
+void BelosBottomSolver<Level>::setOperator(const Teuchos::RCP<matrix_t>& A,
+ Level* level_p)
+{
// make positive definite --> rhs has to change sign as well
A->resumeFill();
@@ -73,7 +77,8 @@ void BelosBottomSolver::setOperator(const Teuchos::RCP<matrix_t>& A) {
}
-std::size_t BelosBottomSolver::getNumIters() {
+template <class Level>
+std::size_t BelosBottomSolver<Level>::getNumIters() {
if ( solver_mp == Teuchos::null )
throw OpalException("BelosBottomSolver::getNumIters()",
"No solver initialized.");
@@ -82,7 +87,8 @@ std::size_t BelosBottomSolver::getNumIters() {
}
-void BelosBottomSolver::initSolver_m(std::string solvertype) {
+template <class Level>
+void BelosBottomSolver<Level>::initSolver_m(std::string solvertype) {
Belos::SolverFactory<scalar_t, mv_t, op_t> factory;
diff --git a/ippl/test/AMR/trilinos/BottomSolver.h b/ippl/test/AMR/trilinos/BottomSolver.h
index b8967438d2f43c9132ee08edc2f861addc9bc9c8..534e198b173b467c346a915a326ef712e33ff6c4 100644
--- a/ippl/test/AMR/trilinos/BottomSolver.h
+++ b/ippl/test/AMR/trilinos/BottomSolver.h
@@ -4,7 +4,7 @@
#include "AmrMultiGridDefs.h"
/// Abstract base class for all base level solvers
-template <class MatrixType, class VectorType>
+template <class Matrix, class Vector, class Level>
class BottomSolver {
public:
@@ -16,14 +16,15 @@ public:
* @param x left-hand side
* @param b right-hand side
*/
- virtual void solve(const VectorType& x,
- const VectorType& b) = 0;
+ virtual void solve(const Vector& x,
+ const Vector& b) = 0;
/*!
* Set the system matrix
* @param A system matrix
*/
- virtual void setOperator(const MatrixType& A) = 0;
+ virtual void setOperator(const Matrix& A,
+ Level* level_p = nullptr) = 0;
/*!
diff --git a/ippl/test/AMR/trilinos/MueLuBottomSolver.h b/ippl/test/AMR/trilinos/MueLuBottomSolver.h
index 6f726813e884afdde8feadc8d3fc8c6781b56842..9790e2c5bdaf928e73d291018e51cd9725cbfb4b 100644
--- a/ippl/test/AMR/trilinos/MueLuBottomSolver.h
+++ b/ippl/test/AMR/trilinos/MueLuBottomSolver.h
@@ -3,14 +3,18 @@
#include "BottomSolver.h"
+#include "Amr/AmrDefs.h"
+
#include <MueLu.hpp>
#include <MueLu_Level.hpp>
#include <MueLu_Utilities.hpp>
#include <MueLu_ParameterListInterpreter.hpp>
+template <class Level>
class MueLuBottomSolver : public BottomSolver<Teuchos::RCP<amr::matrix_t>,
- Teuchos::RCP<amr::multivector_t> >
+ Teuchos::RCP<amr::multivector_t>,
+ Level >
{
public:
typedef amr::matrix_t matrix_t;
@@ -21,6 +25,10 @@ public:
typedef amr::local_ordinal_t lo_t;
typedef amr::global_ordinal_t go_t;
typedef amr::node_t node_t;
+
+ typedef amr::AmrBox_t AmrBox_t;
+ typedef amr::AmrIntVect_t AmrIntVect_t;
+
// typedef amr::AmrGeometry_t AmrGeometry_t;
typedef MueLu::Hierarchy<scalar_t, lo_t, go_t, node_t> hierarchy_t;
@@ -34,12 +42,13 @@ public:
public:
- MueLuBottomSolver(/*const AmrGeometry_t& geom*/);
+ MueLuBottomSolver();
void solve(const Teuchos::RCP<mv_t>& x,
const Teuchos::RCP<mv_t>& b);
- void setOperator(const Teuchos::RCP<matrix_t>& A);
+ void setOperator(const Teuchos::RCP<matrix_t>& A,
+ Level* level_p = nullptr);
std::size_t getNumIters();
@@ -55,10 +64,9 @@ private:
lo_t nSweeps_m; ///< the number of multigrid iterations
- Teuchos::ParameterList mueluList_m;
-
-// const AmrGeometry_t& geom_mr; ///< reference to bottom geometry
-
+ Teuchos::ParameterList mueluList_m;
};
+#include "MueLuBottomSolver.hpp"
+
#endif
diff --git a/ippl/test/AMR/trilinos/MueLuBottomSolver.cpp b/ippl/test/AMR/trilinos/MueLuBottomSolver.hpp
similarity index 57%
rename from ippl/test/AMR/trilinos/MueLuBottomSolver.cpp
rename to ippl/test/AMR/trilinos/MueLuBottomSolver.hpp
index 2e15f4ffc46f533affff4d2b0ab56b84d4280211..185ffdbef7be4689c3de60677258ff22bc08c042 100644
--- a/ippl/test/AMR/trilinos/MueLuBottomSolver.cpp
+++ b/ippl/test/AMR/trilinos/MueLuBottomSolver.hpp
@@ -1,22 +1,23 @@
-#include "MueLuBottomSolver.h"
-
/* https://trilinos.org/wordpress/wp-content/uploads/2015/03/MueLu_tutorial.pdf
* http://prod.sandia.gov/techlib/access-control.cgi/2014/1418624r.pdf
*/
-MueLuBottomSolver::MueLuBottomSolver(/*const AmrGeometry_t& geom*/)
+#include <AMReX.H>
+
+template <class Level>
+MueLuBottomSolver<Level>::MueLuBottomSolver()
: hierarchy_mp(Teuchos::null),
finest_mp(Teuchos::null),
A_mp(Teuchos::null),
- nSweeps_m(6)//,
-// geom_mr(geom)
+ nSweeps_m(4)
{
this->initMueLuList_m();
}
-void MueLuBottomSolver::solve(const Teuchos::RCP<mv_t>& x,
- const Teuchos::RCP<mv_t>& b)
+template <class Level>
+void MueLuBottomSolver<Level>::solve(const Teuchos::RCP<mv_t>& x,
+ const Teuchos::RCP<mv_t>& b)
{
// MueLu requires Xpetra multivectors (wrap them)
Teuchos::RCP<xmv_t> xx = MueLu::TpetraMultiVector_To_XpetraMultiVector(x);
@@ -31,25 +32,45 @@ void MueLuBottomSolver::solve(const Teuchos::RCP<mv_t>& x,
}
-void MueLuBottomSolver::setOperator(const Teuchos::RCP<matrix_t>& A) {
+template <class Level>
+void MueLuBottomSolver<Level>::setOperator(const Teuchos::RCP<matrix_t>& A,
+ Level* level_p) {
A_mp = MueLu::TpetraCrs_To_XpetraMatrix<scalar_t, lo_t, go_t, node_t>(A);
A_mp->SetFixedBlockSize(1); // only 1 DOF per node (pure Laplace problem)
Teuchos::RCP<mv_t> coords_mp = Teuchos::rcp( new amr::multivector_t(A->getDomainMap(),
- 1 /*AMREX_SPACEDIM*/, false) );
-
- const Teuchos::RCP<const amr::dmap_t>& map_r = A->getMap();
-
- // const scalar_t* domain = geom_mr.ProbLo();
- //const scalar_t* spacing = geom_mr.CellSize();
- //for (lo_t d = 0; d < AMREX_SPACEDIM; ++d) {
- for (std::size_t i = 0; i < map_r->getNodeNumElements(); ++i) {
- // coords_mp->replaceLocalValue(i, d, 0.5 * spacing +
- coords_mp->replaceLocalValue(i, 0, map_r->getGlobalElement(i)); //iv[d]);
- }
-// }
-
+ AMREX_SPACEDIM, false) );
+
+ const scalar_t* domain = level_p->geom.ProbLo();
+ const scalar_t* dx = level_p->cellSize();
+#ifdef _OPENMP
+ #pragma omp parallel
+#endif
+ for (amrex::MFIter mfi(level_p->grids, level_p->dmap, true);
+ mfi.isValid(); ++mfi)
+ {
+ const AmrBox_t& tbx = mfi.tilebox();
+ const lo_t* lo = tbx.loVect();
+ const lo_t* hi = tbx.hiVect();
+
+ for (lo_t i = lo[0]; i <= hi[0]; ++i) {
+ for (lo_t j = lo[1]; j <= hi[1]; ++j) {
+#if AMREX_SPACEDIM == 3
+ for (lo_t k = lo[2]; k <= hi[2]; ++k) {
+#endif
+ AmrIntVect_t iv(D_DECL(i, j, k));
+ go_t gidx = level_p->serialize(iv);
+
+ coords_mp->replaceGlobalValue(gidx, 0, (0.5 + i) * dx[0]);
+ coords_mp->replaceGlobalValue(gidx, 1, (0.5 + j) * dx[1]);
+#if AMREX_SPACEDIM == 3
+ coords_mp->replaceGlobalValue(gidx, 2, (0.5 + k) * dx[2]);
+ }
+#endif
+ }
+ }
+ }
Teuchos::RCP<xmv_t> coordinates = MueLu::TpetraMultiVector_To_XpetraMultiVector(coords_mp);
@@ -81,12 +102,14 @@ void MueLuBottomSolver::setOperator(const Teuchos::RCP<matrix_t>& A) {
}
-std::size_t MueLuBottomSolver::getNumIters() {
- return 1;
+template <class Level>
+std::size_t MueLuBottomSolver<Level>::getNumIters() {
+ return nSweeps_m;
}
-void MueLuBottomSolver::initMueLuList_m() {
+template <class Level>
+void MueLuBottomSolver<Level>::initMueLuList_m() {
mueluList_m.set("problem: type", "Poisson-3D");
// mueluList_m.set("problem: symmetric", false);
mueluList_m.set("verbosity", "extreme");
@@ -96,12 +119,15 @@ void MueLuBottomSolver::initMueLuList_m() {
mueluList_m.set("coarse: max size", 200);
mueluList_m.set("multigrid algorithm", "sa");
+ mueluList_m.set("sa: damping factor", 1.33); // default: 1.33
+ mueluList_m.set("sa: use filtered matrix", true);
+ mueluList_m.set("filtered matrix: reuse eigenvalue", false); // false: more expensive
mueluList_m.set("repartition: enable", true);
mueluList_m.set("repartition: rebalance P and R", true);
mueluList_m.set("repartition: partitioner", "zoltan2");
mueluList_m.set("repartition: min rows per proc", 800);
- mueluList_m.set("repartition: start level", 1);
+ mueluList_m.set("repartition: start level", 2);
Teuchos::ParameterList reparms;
reparms.set("algorithm", "rcb");
@@ -111,7 +137,20 @@ void MueLuBottomSolver::initMueLuList_m() {
mueluList_m.set("smoother: type", "CHEBYSHEV");
mueluList_m.set("smoother: pre or post", "both");
-
+ Teuchos::ParameterList smparms;
+ smparms.set("chebyshev: degree", 5);
+ smparms.set("chebyshev: assume matrix does not change", false);
+ smparms.set("chebyshev: zero starting solution", true);
+ mueluList_m.set("smoother: params", smparms);
+
+/*
+ mueluList_m.set("coarse: type", "RELAXATION");
+ Teuchos::ParameterList cparms;
+ cparms.set("relaxation: type", "Gauss-Seidel");
+ cparms.set("relaxation: sweeps", nSweeps_m);
+ cparms.set("relaxation: zero starting solution", true);
+ mueluList_m.set("coarse: params", cparms);
+*/
mueluList_m.set("coarse: type", "KLU2");
mueluList_m.set("aggregation: type", "uncoupled");
@@ -120,28 +159,6 @@ void MueLuBottomSolver::initMueLuList_m() {
mueluList_m.set("transpose: use implicit", false);
- mueluList_m.set("reuse: type", "none");
+ mueluList_m.set("reuse: type", "full"); // none
}
-
-// void MueLuBottomSolver::setupSmoother_m(fManager_t& manager) {
-// std::string ifpackType = "RELAXATION";
-// Teuchos::ParameterList ifpackList;
-// ifpackList.set("relaxation: sweeps", (lo_t) 4);
-// ifpackList.set("relaxation: damping factor", (scalar_t) 1.0);
-// Teuchos::RCP<smootherPrototype_t> smootherPrototype =
-// Teuchos::rcp(new smoother_t(ifpackType, ifpackList));
-//
-// manager.SetFactory("Smoother", Teuchos::rcp(new smootherFactory_t(smootherPrototype)));
-// }
-//
-//
-// void MueLuBottomSolver::setupCoarseSolver_m(fManager_t& manager) {
-// Teuchos::RCP<smootherFactory_t> coarseSolverPrototype =
-// Teuchos::rcp( new solver_t() );
-//
-// Teuchos::RCP<smootherFactory_t> coarseSolverFact =
-// Teuchos::rcp( new smootherFactory_t(coarseSolverPrototype, Teuchos::null) );
-//
-// manager.SetFactory("CoarseSolver", coarseSolverFact);
-// }