MGPoissonSolver.h 8.49 KB
Newer Older
1
//
frey_m's avatar
frey_m committed
2 3 4
// Class MGPoissonSolver
//   This class contains methods for solving Poisson's equation for the
//   space charge portion of the calculation.
5
//
frey_m's avatar
frey_m committed
6 7 8
//   A smoothed aggregation based AMG preconditioned iterative solver for space charge
//   \see FFTPoissonSolver
//   \warning This solver is in an EXPERIMENTAL STAGE. For reliable simulations use the FFTPoissonSolver
9
//
10 11 12
// Copyright (c) 2008,        Yves Ineichen, ETH Zürich,
//               2013 - 2015, Tülin Kaman, Paul Scherrer Institut, Villigen PSI, Switzerland
//               2017 - 2020, Paul Scherrer Institut, Villigen PSI, Switzerland
frey_m's avatar
frey_m committed
13 14
// All rights reserved
//
15 16 17 18 19
// Implemented as part of the master thesis
// "A Parallel Multigrid Solver for Beam Dynamics"
// and the paper
// "A fast parallel Poisson solver on irregular domains applied to beam dynamics simulations"
// (https://doi.org/10.1016/j.jcp.2010.02.022)
frey_m's avatar
frey_m committed
20
//
frey_m's avatar
frey_m committed
21
// In 2020, the code was ported to the second generation Trilinos packages,
frey_m's avatar
frey_m committed
22
// i.e., Epetra --> Tpetra, ML --> MueLu. See also issue #507.
frey_m's avatar
frey_m committed
23
//
frey_m's avatar
frey_m committed
24 25 26 27 28 29 30 31 32
// This file is part of OPAL.
//
// OPAL is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// You should have received a copy of the GNU General Public License
// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
33
//
gsell's avatar
gsell committed
34 35 36 37 38 39 40 41 42
#ifndef MG_POISSON_SOLVER_H_
#define MG_POISSON_SOLVER_H_

//////////////////////////////////////////////////////////////
#include "PoissonSolver.h"
#include "IrregularDomain.h"
//////////////////////////////////////////////////////////////

#include "mpi.h"
43

44 45 46 47 48 49 50 51 52
#include <Tpetra_Vector.hpp>
#include <Tpetra_CrsMatrix.hpp>

#include <Teuchos_DefaultMpiComm.hpp>

#include <BelosLinearProblem.hpp>
#include <BelosSolverFactory.hpp>

#include <MueLu.hpp>
frey_m's avatar
frey_m committed
53 54
#include <MueLu_TpetraOperator.hpp>

55
#include "Teuchos_ParameterList.hpp"
56 57
#include "Algorithms/PartBunch.h"

gsell's avatar
gsell committed
58 59 60 61 62 63 64 65 66 67 68 69 70
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 {
    STD_PREC,
    REUSE_PREC,
    REUSE_HIERARCHY
};

71 72
class BoundaryGeometry;

gsell's avatar
gsell committed
73 74 75
class MGPoissonSolver : public PoissonSolver {

public:
76 77 78 79
    typedef Tpetra::Vector<>                            TpetraVector_t;
    typedef Tpetra::MultiVector<>                       TpetraMultiVector_t;
    typedef Tpetra::Map<>                               TpetraMap_t;
    typedef Tpetra::Vector<>::scalar_type               TpetraScalar_t;
frey_m's avatar
frey_m committed
80
    typedef Tpetra::Vector<>::global_ordinal_type       TpetraGlobalOrdinal_t;
81
    typedef Tpetra::Operator<>                          TpetraOperator_t;
frey_m's avatar
frey_m committed
82
    typedef MueLu::TpetraOperator<>                     MueLuTpetraOperator_t;
83 84 85 86 87 88 89 90 91 92 93 94 95
    typedef Tpetra::CrsMatrix<>                         TpetraCrsMatrix_t;
    typedef Teuchos::MpiComm<int>                       Comm_t;

    typedef Teuchos::ParameterList                      ParameterList_t;

    typedef Belos::SolverManager<TpetraScalar_t,
                                 TpetraMultiVector_t,
                                 TpetraOperator_t>      SolverManager_t;

    typedef Belos::LinearProblem<TpetraScalar_t,
                                 TpetraMultiVector_t,
                                 TpetraOperator_t>      LinearProblem_t;

96 97


frey_m's avatar
frey_m committed
98 99 100 101 102 103
    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);

gsell's avatar
gsell committed
104 105 106 107 108 109 110 111 112 113 114 115
    ~MGPoissonSolver();

    /// given a charge-density field rho and a set of mesh spacings hr, compute
    /// the scalar potential in 'open space'
    /// \param rho (inout) scalar field of the potential
    /// \param hr mesh spacings in each direction
    void computePotential(Field_t &rho, Vector_t hr);
    void computePotential(Field_t &rho, Vector_t hr, double zshift);

    /// set a geometry
    void setGeometry(std::vector<BoundaryGeometry *> geometries);

frey_m's avatar
frey_m committed
116 117 118 119 120 121
    double getXRangeMin(unsigned short /*level*/) { return bp_m->getXRangeMin(); }
    double getXRangeMax(unsigned short /*level*/) { return bp_m->getXRangeMax(); }
    double getYRangeMin(unsigned short /*level*/) { return bp_m->getYRangeMin(); }
    double getYRangeMax(unsigned short /*level*/) { return bp_m->getYRangeMax(); }
    double getZRangeMin(unsigned short /*level*/) { return bp_m->getZRangeMin(); }
    double getZRangeMax(unsigned short /*level*/) { return bp_m->getZRangeMax(); }
122
    void test(PartBunchBase<double, 3>* /*bunch*/) { }
gsell's avatar
gsell committed
123 124 125
    /// useful load balance information
    void printLoadBalanceStats();

126 127
    void extrapolateLHS();

128 129 130 131
    void resizeMesh(Vector_t& origin, Vector_t& hr, const Vector_t& rmin,
                    const Vector_t& rmax, double dh)
    {
        bp_m->resizeMesh(origin, hr, rmin, rmax, dh);
frey_m's avatar
frey_m committed
132 133
    }

gsell's avatar
gsell committed
134 135
    Inform &print(Inform &os) const;

frey_m's avatar
frey_m committed
136 137


gsell's avatar
gsell committed
138 139
private:

frey_m's avatar
fix bug  
frey_m committed
140 141
    bool isMatrixfilled_m;

142 143 144
    // true if CG and GMRES; false if BiCGStab
    bool useLeftPrec_m;

snuverink_j's avatar
snuverink_j committed
145
    //TODO: we need to update this and maybe change attached
gsell's avatar
gsell committed
146 147 148
    //solver!
    /// holding the currently active geometry
    BoundaryGeometry *currentGeometry;
149

gsell's avatar
gsell committed
150 151
    /// container for multiple geometries
    std::vector<BoundaryGeometry *> geometries_m;
152

153
    int repartFreq_m;
gsell's avatar
gsell committed
154 155 156 157 158 159 160 161 162 163 164 165 166 167 168
    /// flag specifying if we are verbose
    bool verbose_m;

    /// tolerance for the iterative solver
    double tol_m;
    /// maximal number of iterations for the iterative solver
    int maxiters_m;
    /// preconditioner mode
    int precmode_m;
    /// maximum number of blocks in Krylov space
    int numBlocks_m;
    /// number of vectors in recycle space
    int recycleBlocks_m;

    /// structure that holds boundary points
frey_m's avatar
frey_m committed
169
    std::unique_ptr<IrregularDomain> bp_m;
gsell's avatar
gsell committed
170 171

    /// right hand side of our problem
172
    Teuchos::RCP<TpetraVector_t> RHS;
gsell's avatar
gsell committed
173
    /// left hand side of the linear system of equations we solve
174
    Teuchos::RCP<TpetraVector_t> LHS;
gsell's avatar
gsell committed
175
    /// matrix used in the linear system of equations
176
    Teuchos::RCP<TpetraCrsMatrix_t> A;
177

frey_m's avatar
frey_m committed
178
    /// Map holding the processor distribution of data
179
    Teuchos::RCP<TpetraMap_t> map_p;
frey_m's avatar
frey_m committed
180

gsell's avatar
gsell committed
181
    /// communicator used by Trilinos
182
    Teuchos::RCP<const Comm_t> comm_mp;
gsell's avatar
gsell committed
183 184

    /// last N LHS's for extrapolating the new LHS as starting vector
gsell's avatar
gsell committed
185
    unsigned int nLHS_m;
186 187
    Teuchos::RCP<TpetraMultiVector_t> P_mp;
    std::deque< TpetraVector_t > OldLHS;
gsell's avatar
gsell committed
188

189 190
    Teuchos::RCP<LinearProblem_t> problem_mp;
    Teuchos::RCP<SolverManager_t>  solver_mp;
191

frey_m's avatar
frey_m committed
192 193 194
    /// MueLu preconditioner object
    Teuchos::RCP<MueLuTpetraOperator_t> prec_mp;

frey_m's avatar
frey_m committed
195
    /// parameter list for the MueLu solver
frey_m's avatar
frey_m committed
196
    Teuchos::ParameterList MueLuList_m;
gsell's avatar
gsell committed
197 198 199
    /// parameter list for the iterative solver (Belos)
    Teuchos::ParameterList belosList;

200 201 202
    /// PartBunch object
    PartBunch *itsBunch_m;

gsell's avatar
gsell committed
203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223
    // mesh and layout objects for rho_m
    Mesh_t *mesh_m;
    FieldLayout_t *layout_m;

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

    /// mesh spacings in each direction
    Vector_t hr_m;
    /// current number of mesh points in each direction
    Vektor<int, 3> nr_m;
    /// global number of mesh points in each direction
    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;
224 225
    IpplTimings::TimerRef FunctionTimer7_m;
    IpplTimings::TimerRef FunctionTimer8_m;
gsell's avatar
gsell committed
226

227 228
    void deletePtr();

frey_m's avatar
frey_m committed
229
    /// recomputes the map
230
    void computeMap(NDIndex<3> localId);
gsell's avatar
gsell committed
231

frey_m's avatar
frey_m committed
232
    /// converts IPPL grid to a 3D map
233 234
    /// \param localId local IPPL grid node indices
    void IPPLToMap3D(NDIndex<3> localId);
gsell's avatar
gsell committed
235 236 237 238 239 240

    /** returns a discretized stencil that has Neumann BC in z direction and
     * Dirichlet BC on the surface of a specified geometry
     * \param hr gridspacings in each direction
     * \param RHS right hand side might be scaled
     */
241
    void ComputeStencil(Vector_t hr, Teuchos::RCP<TpetraVector_t> RHS);
gsell's avatar
gsell committed
242

243 244


gsell's avatar
gsell committed
245 246
protected:
    /// Setup the parameters for the Belos iterative solver.
frey_m's avatar
frey_m committed
247
    void setupBelosList();
gsell's avatar
gsell committed
248 249

    /// Setup the parameters for the SAAMG preconditioner.
frey_m's avatar
frey_m committed
250
    void setupMueLuList();
251 252
};

253

gsell's avatar
gsell committed
254 255 256 257 258
inline Inform &operator<<(Inform &os, const MGPoissonSolver &fs) {
    return fs.print(os);
}

#endif