Branch data 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
|
