MGPoissonSolver.h

Go to the documentation of this file.

// This class contains methods for solving Poisson's equation for the
// space charge portion of the calculation.

#ifdef HAVE_SAAMG_SOLVER

#ifndef MG_POISSON_SOLVER_H_
#define MG_POISSON_SOLVER_H_

#include "PoissonSolver.h"
#include "IrregularDomain.h"
#include "ml_include.h"

#ifdef HAVE_MPI
#include "mpi.h"
#include "Epetra_MpiComm.h"
#else
#include "Epetra_SerialComm.h"
#endif

#if defined(HAVE_ML_EPETRA) && defined(HAVE_ML_TEUCHOS) && defined(HAVE_ML_AZTECOO)

#include "Epetra_Map.h"
#include "Epetra_Vector.h"
#include "Epetra_CrsMatrix.h"
#include "Epetra_MpiComm.h"

#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wunused-local-typedefs"
#include "Teuchos_ParameterList.hpp"
// #include "BelosLinearProblem.hpp"
// #include "BelosRCGSolMgr.hpp"

#include "Algorithms/PartBunch.h"

#include "BelosTypes.hpp"

namespace Belos {
    template <class ScalarType, class MV, class OP>
    class LinearProblem;

    template <class ScalarType, class MV, class OP>
    class OperatorTraits;

    template<class ScalarType, class MV>
    class MultiVecTraits;

    template<class ScalarType, class MV, class OP>
    class SolverManager;

    class EpetraPrecOp;

    template <class ScalarType, class MV, class OP>
    class StatusTestGenResNorm;
}

namespace ML_Epetra {
    class MultiLevelPreconditioner;

    int SetDefaults(std::string ProblemType, Teuchos::ParameterList & List,
                    int * options, double * params, const bool OverWrite);
}

#pragma GCC diagnostic pop

// using Teuchos::RCP;
// using Teuchos::rcp;
// using namespace ML_Epetra;
// using namespace Isorropia;

typedef UniformCartesian<3, double> Mesh_t;
typedef ParticleSpatialLayout<double, 3>::SingleParticlePos_t Vector_t;
typedef ParticleSpatialLayout< double, 3, Mesh_t  > Layout_t;
typedef Cell Center_t;
typedef CenteredFieldLayout<3, Mesh_t, Center_t> FieldLayout_t;
typedef Field<double, 3, Mesh_t, Center_t> Field_t;

enum {
    NO,
    STD_PREC,
    REUSE_PREC,
    REUSE_HIERARCHY
};

class BoundaryGeometry;

class MGPoissonSolver : public PoissonSolver {

public:
    MGPoissonSolver(PartBunch *beam,Mesh_t *mesh, FieldLayout_t *fl, std::vector<BoundaryGeometry *> geometries, std::string itsolver, std::string interpl, double tol, int maxiters, std::string precmode);
    ~MGPoissonSolver();

    void computePotential(Field_t &rho, Vector_t hr);
    void computePotential(Field_t &rho, Vector_t hr, double zshift);

    void setGeometry(std::vector<BoundaryGeometry *> geometries);

    void recomputeMap() { hasParallelDecompositionChanged_m = true; }

    double getXRangeMin(unsigned short level) { return bp->getXRangeMin(); }
    double getXRangeMax(unsigned short level) { return bp->getXRangeMax(); }
    double getYRangeMin(unsigned short level) { return bp->getYRangeMin(); }
    double getYRangeMax(unsigned short level) { return bp->getYRangeMax(); }
    double getZRangeMin(unsigned short level) { return bp->getZRangeMin(); }
    double getZRangeMax(unsigned short level) { return bp->getZRangeMax(); }
    void test(PartBunchBase<double, 3> *bunch) { }
    void printLoadBalanceStats();

    void extrapolateLHS();

    Inform &print(Inform &os) const;

private:

    //TODO: we need to update this and maybe change attached
    //solver!
    BoundaryGeometry *currentGeometry;

    std::vector<BoundaryGeometry *> geometries_m;

    bool hasGeometryChanged_m;
    bool hasParallelDecompositionChanged_m;
    int repartFreq_m;
    bool useRCB_m;
    bool verbose_m;

    double tol_m;
    int maxiters_m;
    int itsolver_m;
    int precmode_m;
    int numIter_m;
    int tolerableIterationsCount_m;
    bool forcePreconditionerRecomputation_m;
    int numBlocks_m;
    int recycleBlocks_m;

    IrregularDomain *bp;

    Teuchos::RCP<Epetra_Vector> RHS;
    Teuchos::RCP<Epetra_Vector> LHS;
    Teuchos::RCP<Epetra_CrsMatrix> A;

    ML_Epetra::MultiLevelPreconditioner *MLPrec;
    Epetra_Map *Map;
    Epetra_MpiComm Comm;

    uint nLHS_m;
    Teuchos::RCP<Epetra_MultiVector> P;
    std::deque< Epetra_Vector > OldLHS;

    typedef double                          ST;
    typedef Epetra_Operator                 OP;
    typedef Epetra_MultiVector              MV;
    typedef Belos::OperatorTraits<ST, MV, OP> OPT;
    typedef Belos::MultiVecTraits<ST, MV>    MVT;

    //Belos::LinearProblem<double, MV, OP> problem;
    typedef Belos::LinearProblem<ST, MV, OP> problem;
    Teuchos::RCP< problem > problem_ptr;

    typedef Belos::SolverManager<ST, MV, OP> solver;
    Teuchos::RCP< solver > solver_ptr;

    Teuchos::RCP< Belos::EpetraPrecOp > prec_m;
    Teuchos::RCP< Belos::StatusTestGenResNorm< ST, MV, OP > > convStatusTest;

    Teuchos::ParameterList MLList_m;
    Teuchos::ParameterList belosList;

    PartBunch *itsBunch_m;

    // mesh and layout objects for rho_m
    Mesh_t *mesh_m;
    FieldLayout_t *layout_m;

    // domains for the various fields
    NDIndex<3> domain_m;

    Vector_t hr_m;
    Vektor<int, 3> nr_m;
    Vektor<int, 3> orig_nr_m;

    // timers
    IpplTimings::TimerRef FunctionTimer1_m;
    IpplTimings::TimerRef FunctionTimer2_m;
    IpplTimings::TimerRef FunctionTimer3_m;
    IpplTimings::TimerRef FunctionTimer4_m;
    IpplTimings::TimerRef FunctionTimer5_m;
    IpplTimings::TimerRef FunctionTimer6_m;
    IpplTimings::TimerRef FunctionTimer7_m;
    IpplTimings::TimerRef FunctionTimer8_m;

    void deletePtr();

    void computeMap(NDIndex<3> localId);

    void redistributeWithRCB(NDIndex<3> localId);

    void IPPLToMap3D(NDIndex<3> localId);

    void ComputeStencil(Vector_t hr, Teuchos::RCP<Epetra_Vector> RHS);



protected:

    inline void SetupBelosList() {
        belosList.set("Maximum Iterations", maxiters_m);
        belosList.set("Convergence Tolerance", tol_m);

        if(numBlocks_m != 0 && recycleBlocks_m != 0){//only set if solver==RCGSolMgr
            belosList.set("Num Blocks", numBlocks_m);               // Maximum number of blocks in Krylov space
            belosList.set("Num Recycled Blocks", recycleBlocks_m); // Number of vectors in recycle space
        }
        if(verbose_m) {
            belosList.set("Verbosity", Belos::Errors + Belos::Warnings + Belos::TimingDetails + Belos::FinalSummary + Belos::StatusTestDetails);
            belosList.set("Output Frequency", 1);
        } else
            belosList.set("Verbosity", Belos::Errors + Belos::Warnings);
      }

    inline void SetupMLList() {
        ML_Epetra::SetDefaults("SA", MLList_m, 0, 0, true);

        MLList_m.set("max levels", 8);
        MLList_m.set("increasing or decreasing", "increasing");

        // we use a V-cycle
        MLList_m.set("prec type", "MGV");

        // uncoupled aggregation is used (every processor aggregates
        // only local data)
        MLList_m.set("aggregation: type", "Uncoupled");

        // smoother related parameters
        MLList_m.set("smoother: type","Chebyshev");
        MLList_m.set("smoother: sweeps", 3);
        MLList_m.set("smoother: pre or post", "both");

        // on the coarsest level we solve with  Tim Davis' implementation of
        // Gilbert-Peierl's left-looking sparse partial pivoting algorithm,
        // with Eisenstat & Liu's symmetric pruning. Gilbert's version appears
        // as \c [L,U,P]=lu(A) in MATLAB. It doesn't exploit dense matrix
        // kernels, but it is the only sparse LU factorization algorithm known to be
        // asymptotically optimal, in the sense that it takes time proportional to the
        // number of floating-point operations.
        MLList_m.set("coarse: type", "Amesos-KLU");

        //XXX: or use Chebyshev coarse level solver
        // SEE PAPER FOR EVALUATION KLU vs. Chebyshev
        //MLList.set("coarse: sweeps", 10);
        //MLList.set("coarse: type", "Chebyshev");

        // Controls the amount of printed information.
        // Ranges from 0 to 10 (0 is no output, and
        // 10 is incredibly detailed output). Default: 0
        if(verbose_m)
            MLList_m.set("ML output", 10);

        // heuristic for max coarse size depending on number of processors
        int coarsest_size = std::max(Comm.NumProc() * 10, 1024);
        MLList_m.set("coarse: max size", coarsest_size);

    }

};



inline Inform &operator<<(Inform &os, const MGPoissonSolver &fs) {
    return fs.print(os);
}

#else

#include <stdlib.h>
#include <stdio.h>
#ifdef HAVE_MPI
#include "mpi.h"
#endif

int main(int argc, char *argv[]) {
#ifdef HAVE_MPI
    MPI_Init(&argc, &argv);
#endif

    puts("Please configure ML with:");
    puts("--enable-epetra");
    puts("--enable-teuchos");
    puts("--enable-aztecoo");

#ifdef HAVE_MPI
    MPI_Finalize();
#endif

    return(EXIT_SUCCESS);
}

#endif /* #if defined(HAVE_ML_EPETRA) && defined(HAVE_ML_TEUCHOS) && defined(HAVE_ML_AZTECOO) */

#endif /* #ifndef MG_POISSON_SOLVER_H_ */

#endif /* #ifdef HAVE_SAAMG_SOLVER */