AMR-Trilinos: Preconditioner update

ippl/test/AMR/trilinos/AmrPreconditioner.h

@@ -12,7 +12,7 @@ namespace amr {
         SA,             ///< smoothed aggregation multigrid
         JACOBI,         ///< Jacobi point relaxation
         BLOCK_JACOBI,   ///< Jacobi block relaxation
-        GS,             ///> Gauss-Seidel point relaxation
+        GS,             ///< Gauss-Seidel point relaxation
         BLOCK_GS        ///< Gauss-Seidel block relaxation
     };
 }

src/Solvers/AMR_MG/AmrPreconditioner.h

@@ -12,13 +12,13 @@ namespace amr {
         SA,             ///< smoothed aggregation multigrid
         JACOBI,         ///< Jacobi point relaxation
         BLOCK_JACOBI,   ///< Jacobi block relaxation
-        GS,             ///> Gauss-Seidel point relaxation
+        GS,             ///< Gauss-Seidel point relaxation
         BLOCK_GS        ///< Gauss-Seidel block relaxation
     };
 }
 
 /// Bottom solver preconditioners
-template <class MatrixType>
+template <class Matrix, class Level>
 class AmrPreconditioner
 {
 public:
@@ -29,9 +29,13 @@ public:
     /*!
      * Instantiate the preconditioner matrix
      * @param A matrix for which to create preconditioner
+     * @param level_p bottom level if necessary to build preconditioner
      */
-    virtual void create(const Teuchos::RCP<MatrixType>& A) = 0;
+    virtual void create(const Teuchos::RCP<Matrix>& A, Level* level_p = nullptr) = 0;
     
+    /*!
+     * @returns the preconditioner
+     */
     virtual Teuchos::RCP<operator_t> get() = 0;
 };
 

src/Solvers/AMR_MG/CMakeLists.txt

 set (_SRCS
-    AmesosBottomSolver.cpp
     AmrMultiGrid.cpp
     AmrSmoother.cpp
-    BelosBottomSolver.cpp
-    Ifpack2Preconditioner.cpp
-    MueLuPreconditioner.cpp
 )
 
 include_directories (
@@ -15,6 +11,7 @@ add_opal_sources(${_SRCS})
 
 set (HDRS
     AmesosBottomSolver.h
+    AmesosBottomSolver.hpp
     AmrBoundary.h
     AmrDirichletBoundary.h
     AmrInterpolater.h
@@ -34,9 +31,13 @@ set (HDRS
     AmrTrilinearInterpolater.h
     AmrTrilinearInterpolater.hpp
     BelosBottomSolver.h
+    BelosBottomSolver.hpp
     BottomSolver.h
     Ifpack2Preconditioner.h
+    Ifpack2Preconditioner.hpp
     MueLuPreconditioner.h
+    MueLuPreconditioner.hpp
 )
 
-install (FILES ${HDRS} DESTINATION "${CMAKE_INSTALL_PREFIX}/include/Solvers/AMR_MG/")
\ No newline at end of file
+
+install (FILES ${HDRS} DESTINATION "${CMAKE_INSTALL_PREFIX}/include/Solvers/AMR_MG/")

src/Solvers/AMR_MG/Ifpack2Preconditioner.h

@@ -5,7 +5,8 @@
 
 #include <Ifpack2_Factory.hpp>
 
-class Ifpack2Preconditioner : public AmrPreconditioner<amr::matrix_t>
+template <class Level>
+class Ifpack2Preconditioner : public AmrPreconditioner<amr::matrix_t, Level>
 {
 public:
     typedef amr::Preconditioner Preconditioner;
@@ -22,7 +23,7 @@ public:
     
     Ifpack2Preconditioner(Preconditioner prec);
     
-    void create(const Teuchos::RCP<amr::matrix_t>& A);
+    void create(const Teuchos::RCP<amr::matrix_t>& A, Level* level_p = nullptr);
     
     Teuchos::RCP<amr::operator_t> get();
     
@@ -44,4 +45,6 @@ private:
     Teuchos::RCP<precond_t> prec_mp;
 };
 
+#include "Ifpack2Preconditioner.hpp"
+
 #endif

src/Solvers/AMR_MG/Ifpack2Preconditioner.cpp → src/Solvers/AMR_MG/Ifpack2Preconditioner.hpp

-#include "Ifpack2Preconditioner.h"
-
 #include "Utilities/OpalException.h"
 
-Ifpack2Preconditioner::Ifpack2Preconditioner(Preconditioner prec)
+template <class Level>
+Ifpack2Preconditioner<Level>::Ifpack2Preconditioner(Preconditioner prec)
     : prec_mp(Teuchos::null)
 {
     this->init_m(prec);
 }
 
 
-void Ifpack2Preconditioner::create(const Teuchos::RCP<amr::matrix_t>& A) {
+template <class Level>
+void Ifpack2Preconditioner<Level>::create(const Teuchos::RCP<amr::matrix_t>& A,
+                                          Level* level_p)
+{
     Ifpack2::Factory factory;
     
     prec_mp = factory.create(prectype_m, A.getConst());
@@ -19,12 +21,14 @@ void Ifpack2Preconditioner::create(const Teuchos::RCP<amr::matrix_t>& A) {
 }
 
 
-Teuchos::RCP<amr::operator_t> Ifpack2Preconditioner::get() {
+template <class Level>
+Teuchos::RCP<amr::operator_t> Ifpack2Preconditioner<Level>::get() {
     return prec_mp;
 }
 
 
-void Ifpack2Preconditioner::fillMap(map_t& map) {
+template <class Level>
+void Ifpack2Preconditioner<Level>::fillMap(map_t& map) {
     map["ILUT"]         = Preconditioner::ILUT;
     map["CHEBYSHEV"]    = Preconditioner::CHEBYSHEV;
     map["RILUK"]        = Preconditioner::RILUK;
@@ -35,7 +39,8 @@ void Ifpack2Preconditioner::fillMap(map_t& map) {
 }
 
 
-void Ifpack2Preconditioner::init_m(Preconditioner prec)
+template <class Level>
+void Ifpack2Preconditioner<Level>::init_m(Preconditioner prec)
 {
     params_mp = Teuchos::parameterList();
     

src/Solvers/AMR_MG/MueLuPreconditioner.h

@@ -7,15 +7,21 @@
 #include <MueLu.hpp>
 #include <MueLu_TpetraOperator.hpp>
 
-class MueLuPreconditioner : public AmrPreconditioner<amr::matrix_t>
+template <class Level>
+class MueLuPreconditioner : public AmrPreconditioner<amr::matrix_t, Level>
 {
 public:
     typedef amr::Preconditioner Preconditioner;
     
+    typedef amr::scalar_t scalar_t;
+    typedef amr::local_ordinal_t lo_t;
+    typedef amr::global_ordinal_t go_t;
+    typedef amr::AmrBox_t AmrBox_t;
+    
     typedef MueLu::TpetraOperator<
-        amr::scalar_t,
-        amr::local_ordinal_t,
-        amr::global_ordinal_t,
+        scalar_t,
+        lo_t,
+        go_t,
         amr::node_t
     > precond_t;
     
@@ -25,10 +31,9 @@ public:
     
 public:
     
-    MueLuPreconditioner(const AmrIntVect_t& grid,
-                        const bool& rebalance);
+    MueLuPreconditioner(const bool& rebalance);
     
-    void create(const Teuchos::RCP<amr::matrix_t>& A);
+    void create(const Teuchos::RCP<amr::matrix_t>& A, Level* level_p =  nullptr);
     
     Teuchos::RCP<amr::operator_t> get();
     
@@ -43,9 +48,10 @@ private:
     Teuchos::RCP<precond_t> prec_mp;
 
     Teuchos::RCP<amr::multivector_t> coords_mp;
-
-    const AmrIntVect_t grid_m;
+    
     const bool rebalance_m;
 };
 
+#include "MueLuPreconditioner.hpp"
+
 #endif

src/Solvers/AMR_MG/MueLuPreconditioner.cpp → src/Solvers/AMR_MG/MueLuPreconditioner.hpp

-#include "MueLuPreconditioner.h"
+#include <AMReX.H>
 
 #include <MueLu_CreateTpetraPreconditioner.hpp>
 
-
-
-MueLuPreconditioner::MueLuPreconditioner(const AmrIntVect_t& grid,
-                                         const bool& rebalance)
+template <class Level>
+MueLuPreconditioner<Level>::MueLuPreconditioner(const bool& rebalance)
     : prec_mp(Teuchos::null),
       coords_mp(Teuchos::null),
-      grid_m(grid),
       rebalance_m(rebalance)
 {
     this->init_m();
 }
 
 
-void MueLuPreconditioner::create(const Teuchos::RCP<amr::matrix_t>& A) {
+template <class Level>
+void MueLuPreconditioner<Level>::create(const Teuchos::RCP<amr::matrix_t>& A,
+                                        Level* level_p)
+{
 
     coords_mp = Teuchos::null;
 
     if ( rebalance_m ) {
-        coords_mp = Teuchos::rcp( new amr::multivector_t(A->getDomainMap(),
-                                                         1, false) );
-    
-        const Teuchos::RCP<const amr::dmap_t>& map_r = A->getMap();
+        
+        const scalar_t* domain = level_p->geom.ProbLo();
+        const scalar_t* dx = level_p->cellSize();
+        
+        coords_mp = Teuchos::rcp(
+            new amr::multivector_t(A->getDomainMap(), AMREX_SPACEDIM, false)
+        );
     
-        for (std::size_t i = 0; i < map_r->getNodeNumElements(); ++i)
-            coords_mp->replaceLocalValue(i, 0, map_r->getGlobalElement(i));
+#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, domain[0] + (0.5 + i) * dx[0]);
+                        coords_mp->replaceGlobalValue(gidx, 1, domain[1] + (0.5 + j) * dx[1]);
+#if AMREX_SPACEDIM == 3
+                        coords_mp->replaceGlobalValue(gidx, 2, domain[2] + (0.5 + k) * dx[2]);
+                    }
+#endif
+                }
+            }
+        }
     }
 
     prec_mp = MueLu::CreateTpetraPreconditioner(A, params_m, coords_mp);
@@ -35,19 +63,21 @@ void MueLuPreconditioner::create(const Teuchos::RCP<amr::matrix_t>& A) {
 }
 
 
-Teuchos::RCP<amr::operator_t> MueLuPreconditioner::get() {
+template <class Level>
+Teuchos::RCP<amr::operator_t> MueLuPreconditioner<Level>::get() {
     return prec_mp;
 }
 
 
-void MueLuPreconditioner::fillMap(map_t& map) {
+template <class Level>
+void MueLuPreconditioner<Level>::fillMap(map_t& map) {
     map["SA"] = Preconditioner::SA;
 }
 
 
-void MueLuPreconditioner::init_m() {
+template <class Level>
+void MueLuPreconditioner<Level>::init_m() {
     params_m.set("problem: type", "Poisson-3D");
-//    params_m.set("problem: symmetric", false);
     params_m.set("verbosity", "extreme");
     params_m.set("number of equations", 1);
     params_m.set("max levels", 8);
@@ -55,6 +85,9 @@ void MueLuPreconditioner::init_m() {
 
     params_m.set("coarse: max size", 200);
     params_m.set("multigrid algorithm", "sa");
+    params_m.set("sa: damping factor", 1.33); // default: 1.33
+    params_m.set("sa: use filtered matrix", true);
+    params_m.set("filtered matrix: reuse eigenvalue", false); // false: more expensive
     
     params_m.set("repartition: enable", rebalance_m);
     params_m.set("repartition: rebalance P and R", rebalance_m);
@@ -63,11 +96,20 @@ void MueLuPreconditioner::init_m() {
     params_m.set("repartition: start level", 2);
 
     Teuchos::ParameterList reparms;
-    reparms.set("algorithm", "rcb");
+    reparms.set("algorithm", "multijagged"); // rcb
     //    reparms.set("partitioning_approach", "partition");
 
     params_m.set("repartition: params", reparms);
     
+    
+    params_m.set("smoother: type", "CHEBYSHEV");
+    params_m.set("smoother: pre or post", "both");
+    Teuchos::ParameterList smparms;
+    smparms.set("chebyshev: degree", 3);
+    smparms.set("chebyshev: assume matrix does not change", false);
+    smparms.set("chebyshev: zero starting solution", true);
+    params_m.set("smoother: params", smparms);
+    
     params_m.set("smoother: type", "CHEBYSHEV");
     params_m.set("smoother: pre or post", "both");
 
@@ -77,7 +119,7 @@ void MueLuPreconditioner::init_m() {
     params_m.set("aggregation: min agg size", 3);
     params_m.set("aggregation: max agg size", 27);
 
-    params_m.set("transpose: use implicit", true);
+    params_m.set("transpose: use implicit", false);
 
-    params_m.set("reuse: type", "none");
+    params_m.set("reuse: type", "full"); // none
 }