AMR-Trilinos: add new files in OPAL

src/Solvers/AMR_MG/CMakeLists.txt

@@ -35,6 +35,8 @@ set (HDRS
     BottomSolver.h
     Ifpack2Preconditioner.h
     Ifpack2Preconditioner.hpp
+    MueLuBottomSolver.h
+    MueLuBottomSolver.hpp
     MueLuPreconditioner.h
     MueLuPreconditioner.hpp
 )

src/Solvers/AMR_MG/MueLuBottomSolver.h

+#ifndef MUELU_BOTTOM_SOLVER_H
+#define MUELU_BOTTOM_SOLVER_H
+
+#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>,
+                                              Level >
+{
+public:
+    typedef amr::matrix_t matrix_t;
+    typedef amr::vector_t vector_t;
+    typedef amr::scalar_t scalar_t;
+    typedef amr::multivector_t mv_t;
+    typedef amr::operator_t op_t;
+    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;
+    typedef MueLu::Level level_t;
+    typedef Xpetra::Matrix<scalar_t, lo_t, go_t, node_t> xmatrix_t;
+    typedef Xpetra::MultiVector<scalar_t, lo_t, go_t, node_t> xmv_t;
+    typedef MueLu::Utilities<scalar_t, lo_t, go_t, node_t> util_t;
+    
+    typedef MueLu::ParameterListInterpreter<scalar_t, lo_t, go_t, node_t> pListInterpreter_t;
+    typedef MueLu::HierarchyManager<scalar_t, lo_t, go_t, node_t> manager_t;
+        
+public:
+    
+    MueLuBottomSolver();
+    
+    void solve(const Teuchos::RCP<mv_t>& x,
+               const Teuchos::RCP<mv_t>& b);
+    
+    void setOperator(const Teuchos::RCP<matrix_t>& A,
+                     Level* level_p = nullptr);
+    
+    std::size_t getNumIters();
+
+private:
+    void initMueLuList_m();
+    
+private:
+    Teuchos::RCP<hierarchy_t> hierarchy_mp;     ///< manages the multigrid hierarchy
+    
+    Teuchos::RCP<level_t> finest_mp;            ///< finest level of hierarchy
+
+    Teuchos::RCP<xmatrix_t> A_mp;               ///< MueLu requires Xpetra
+
+    lo_t nSweeps_m;                             ///< the number of multigrid iterations
+    
+    Teuchos::ParameterList mueluList_m;    
+};
+
+#include "MueLuBottomSolver.hpp"
+
+#endif

src/Solvers/AMR_MG/MueLuBottomSolver.hpp

+/* https://trilinos.org/wordpress/wp-content/uploads/2015/03/MueLu_tutorial.pdf
+ * http://prod.sandia.gov/techlib/access-control.cgi/2014/1418624r.pdf
+ */
+
+#include <AMReX.H>
+
+template <class Level>
+MueLuBottomSolver<Level>::MueLuBottomSolver()
+    : hierarchy_mp(Teuchos::null),
+      finest_mp(Teuchos::null),
+      A_mp(Teuchos::null),
+      nSweeps_m(4)
+{
+    this->initMueLuList_m();
+}
+
+
+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);
+    Teuchos::RCP<xmv_t> xb = MueLu::TpetraMultiVector_To_XpetraMultiVector(b);
+
+    // InitialGuessIsZero = true
+    // startLevel = 0
+    hierarchy_mp->Iterate(*xb, *xx, nSweeps_m, true, 0);
+    
+    // put multivector back
+    x->assign(*util_t::MV2NonConstTpetraMV2(*xx));
+}
+
+
+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(), 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, 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
+            }
+        }
+    }
+
+    Teuchos::RCP<xmv_t> coordinates = MueLu::TpetraMultiVector_To_XpetraMultiVector(coords_mp);
+
+
+    Teuchos::RCP<manager_t> mueluFactory = Teuchos::rcp(
+        new pListInterpreter_t(mueluList_m)
+    );
+
+    // empty multigrid hierarchy with a finest level only
+    hierarchy_mp = mueluFactory->CreateHierarchy();
+
+    hierarchy_mp->GetLevel(0)->Set("A", A_mp);
+
+    Teuchos::RCP<mv_t> nullspace = Teuchos::rcp(new mv_t(A->getRowMap(), 1));
+    Teuchos::RCP<xmv_t> xnullspace = MueLu::TpetraMultiVector_To_XpetraMultiVector(nullspace);
+    xnullspace->putScalar(1.0);
+
+    hierarchy_mp->GetLevel(0)->Set("Nullspace", xnullspace);
+    hierarchy_mp->GetLevel(0)->Set("Coordinates", coordinates);
+    hierarchy_mp->IsPreconditioner(false);
+    hierarchy_mp->setDefaultVerbLevel(Teuchos::VERB_HIGH);
+   
+    finest_mp = hierarchy_mp->GetLevel();
+    finest_mp->Set("A", A_mp);
+    finest_mp->Set("Coordinates", coordinates);
+    finest_mp->Set("Nullspace", xnullspace);
+    
+    mueluFactory->SetupHierarchy(*hierarchy_mp);
+}
+
+
+template <class Level>
+std::size_t MueLuBottomSolver<Level>::getNumIters() {
+    return nSweeps_m;
+}
+
+
+template <class Level>
+void MueLuBottomSolver<Level>::initMueLuList_m() {
+    mueluList_m.set("problem: type", "Poisson-3D");
+    mueluList_m.set("verbosity", "extreme");
+    mueluList_m.set("number of equations", 1);
+    mueluList_m.set("max levels", 8);
+    mueluList_m.set("cycle type", "V");
+
+    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", 2);
+
+    Teuchos::ParameterList reparms;
+    reparms.set("algorithm", "multijagged"); // rcb
+    //    reparms.set("partitioning_approach", "partition");
+
+    mueluList_m.set("repartition: params", reparms);
+    
+    mueluList_m.set("smoother: type", "CHEBYSHEV");
+    mueluList_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);
+    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);
+    cparms.set("relaxation: use l1", true);
+    cparms.set("relaxation: l1 eta", 1.5);
+    mueluList_m.set("coarse: params", cparms);
+
+//    mueluList_m.set("coarse: type", "KLU2");
+
+    mueluList_m.set("aggregation: type", "uncoupled");
+    mueluList_m.set("aggregation: min agg size", 3);
+    mueluList_m.set("aggregation: max agg size", 27);
+
+    mueluList_m.set("transpose: use implicit", false);
+
+    mueluList_m.set("reuse: type", "full"); // none
+}
+