Line data    Source code

       1              : //
       2              : // Class P3MSolver
       3              : //   Poisson solver for periodic boundaries, based on FFTs.
       4              : //   Solves laplace(phi) = -rho, and E = -grad(phi).
       5              : //
       6              : //   Uses a convolution with a Green's function given by:
       7              : //      G(r) = ke * erf(alpha * r) / r,
       8              : //   where ke = Coulomb constant,
       9              : //         alpha = controls long-range interaction.
      10              : //
      11              : //
      12              : 
      13              : #ifndef IPPL_P3M_SOLVER_H_
      14              : #define IPPL_P3M_SOLVER_H_
      15              : 
      16              : #include "Types/Vector.h"
      17              : 
      18              : #include "Field/Field.h"
      19              : 
      20              : #include "FFT/FFT.h"
      21              : #include "FieldLayout/FieldLayout.h"
      22              : #include "Meshes/UniformCartesian.h"
      23              : #include "Poisson.h"
      24              : 
      25              : namespace ippl {
      26              :     template <typename FieldLHS, typename FieldRHS>
      27              :     class P3MSolver : public Poisson<FieldLHS, FieldRHS> {
      28              :         constexpr static unsigned Dim = FieldLHS::dim;
      29              :         using Trhs                    = typename FieldRHS::value_type;
      30              :         using mesh_type               = typename FieldRHS::Mesh_t;
      31              : 
      32              :     public:
      33              :         // type of output
      34              :         using Base = Poisson<FieldLHS, FieldRHS>;
      35              : 
      36              :         // types for LHS and RHS
      37              :         using typename Base::lhs_type, typename Base::rhs_type;
      38              : 
      39              :         // define a type for the 3 dimensional real to complex Fourier transform
      40              :         typedef FFT<RCTransform, FieldRHS> FFT_t;
      41              : 
      42              :         // define a type for a 3 dimensional field (e.g. charge density field)
      43              :         // define a type of Field with integers to be used for the helper Green's function
      44              :         // also define a type for the Fourier transformed complex valued fields
      45              :         typedef FieldRHS Field_t;
      46              :         typedef Field<int, Dim, mesh_type, typename FieldLHS::Centering_t> IField_t;
      47              :         typedef typename FFT_t::ComplexField CxField_t;
      48              :         typedef Vector<Trhs, Dim> Vector_t;
      49              : 
      50              :         // define type for field layout
      51              :         typedef FieldLayout<Dim> FieldLayout_t;
      52              : 
      53              :         // constructor and destructor
      54              :         P3MSolver();
      55              :         P3MSolver(rhs_type& rhs, ParameterList& params);
      56              :         P3MSolver(lhs_type& lhs, rhs_type& rhs, ParameterList& params);
      57            0 :         ~P3MSolver() = default;
      58              : 
      59              :         // override the setRhs function of the Solver class
      60              :         // since we need to call initializeFields()
      61              :         void setRhs(rhs_type& rhs) override;
      62              : 
      63              :         // solve the Poisson equation
      64              :         // more specifically, compute the scalar potential given a density field rho
      65              :         void solve() override;
      66              : 
      67              :         // function called in the constructor to initialize the fields
      68              :         void initializeFields();
      69              : 
      70              :         // compute standard Green's function
      71              :         void greensFunction();
      72              : 
      73              :     private:
      74              :         Field_t grn_m;  // the Green's function
      75              : 
      76              :         CxField_t rhotr_m;
      77              :         CxField_t grntr_m;
      78              :         CxField_t tempFieldComplex_m;
      79              : 
      80              :         // fields that facilitate the calculation in greensFunction()
      81              :         IField_t grnIField_m[Dim];
      82              : 
      83              :         // the FFT object
      84              :         std::unique_ptr<FFT_t> fft_m;
      85              : 
      86              :         // mesh and layout objects for rho_m (RHS)
      87              :         mesh_type* mesh_mp;
      88              :         FieldLayout_t* layout_mp;
      89              : 
      90              :         // mesh and layout objects for the Fourier transformed Complex fields
      91              :         std::unique_ptr<mesh_type> meshComplex_m;
      92              :         std::unique_ptr<FieldLayout_t> layoutComplex_m;
      93              : 
      94              :         // domains for the various fields
      95              :         NDIndex<Dim> domain_m;         // physical domain
      96              :         NDIndex<Dim> domainComplex_m;  // Fourier domain
      97              : 
      98              :         // mesh spacing and mesh size
      99              :         Vector_t hr_m;
     100              :         Vector<int, Dim> nr_m;
     101              : 
     102              :     protected:
     103            0 :         virtual void setDefaultParameters() override {
     104              :             using heffteBackend       = typename FFT_t::heffteBackend;
     105            0 :             heffte::plan_options opts = heffte::default_options<heffteBackend>();
     106            0 :             this->params_m.add("use_pencils", opts.use_pencils);
     107            0 :             this->params_m.add("use_reorder", opts.use_reorder);
     108            0 :             this->params_m.add("use_gpu_aware", opts.use_gpu_aware);
     109            0 :             this->params_m.add("r2c_direction", 0);
     110              : 
     111            0 :             switch (opts.algorithm) {
     112            0 :                 case heffte::reshape_algorithm::alltoall:
     113            0 :                     this->params_m.add("comm", a2a);
     114            0 :                     break;
     115            0 :                 case heffte::reshape_algorithm::alltoallv:
     116            0 :                     this->params_m.add("comm", a2av);
     117            0 :                     break;
     118            0 :                 case heffte::reshape_algorithm::p2p:
     119            0 :                     this->params_m.add("comm", p2p);
     120            0 :                     break;
     121            0 :                 case heffte::reshape_algorithm::p2p_plined:
     122            0 :                     this->params_m.add("comm", p2p_pl);
     123            0 :                     break;
     124            0 :                 default:
     125            0 :                     throw IpplException("P3MSolver::setDefaultParameters",
     126              :                                         "Unrecognized heffte communication type");
     127              :             }
     128            0 :         }
     129              :     };
     130              : }  // namespace ippl
     131              : 
     132              : #include "PoissonSolvers/P3MSolver.hpp"
     133              : #endif