Branch data Line data Source code
1 : : //
2 : : // Helper functions used in PoissonCG.h
3 : : //
4 : :
5 : : #ifndef IPPL_LAPLACE_HELPERS_H
6 : : #define IPPL_LAPLACE_HELPERS_H
7 : : namespace ippl {
8 : : namespace detail {
9 : : // Implements the poisson matrix acting on a d dimensional field
10 : : template <typename E>
11 : : struct meta_poisson : public Expression<meta_poisson<E>, sizeof(E)> {
12 : : constexpr static unsigned dim = E::dim;
13 : :
14 : : KOKKOS_FUNCTION
15 : : meta_poisson(const E& u)
16 : : : u_m(u) {}
17 : :
18 : : template <typename... Idx>
19 : : KOKKOS_INLINE_FUNCTION auto operator()(const Idx... args) const {
20 : : using index_type = std::tuple_element_t<0, std::tuple<Idx...>>;
21 : : using T = typename E::Mesh_t::value_type;
22 : : T res = 0;
23 : : index_type coords[dim] = {args...};
24 : : auto&& center = apply(u_m, coords);
25 : : for (unsigned d = 0; d < dim; d++) {
26 : : index_type coords[dim] = {args...};
27 : :
28 : : coords[d] -= 1;
29 : : auto&& left = apply(u_m, coords);
30 : :
31 : : coords[d] += 2;
32 : : auto&& right = apply(u_m, coords);
33 : : res += (2.0 * center - left - right);
34 : : }
35 : : return res;
36 : : }
37 : :
38 : : private:
39 : : const E u_m;
40 : : };
41 : :
42 : : template <typename E>
43 : : struct meta_lower_laplace
44 : : : public Expression<meta_lower_laplace<E>,
45 : : sizeof(E) + sizeof(typename E::Mesh_t::vector_type)
46 : : + 2 * sizeof(typename E::Layout_t::NDIndex_t)
47 : : + sizeof(unsigned)> {
48 : : constexpr static unsigned dim = E::dim;
49 : :
50 : : KOKKOS_FUNCTION
51 : 12 : meta_lower_laplace(const E& u, const typename E::Mesh_t::vector_type& hvector,
52 : : unsigned nghosts, const typename E::Layout_t::NDIndex_t& ldom,
53 : : const typename E::Layout_t::NDIndex_t& domain)
54 : 12 : : u_m(u)
55 : 12 : , hvector_m(hvector)
56 : 12 : , nghosts_m(nghosts)
57 : 12 : , ldom_m(ldom)
58 : 12 : , domain_m(domain) {}
59 : :
60 : : /*
61 : : * n-dimensional lower triangular Laplacian
62 : : */
63 : : template <typename... Idx>
64 : 238656 : KOKKOS_INLINE_FUNCTION auto operator()(const Idx... args) const {
65 : : using index_type = std::tuple_element_t<0, std::tuple<Idx...>>;
66 : : using T = typename E::Mesh_t::value_type;
67 : 238656 : T res = 0;
68 : :
69 [ + + ]: 1510528 : for (unsigned d = 0; d < dim; d++) {
70 : 1271872 : index_type coords[dim] = {args...};
71 : 1271872 : const int global_index = coords[d] + ldom_m[d].first() - nghosts_m;
72 [ + - ]: 1271872 : const int size = domain_m.length()[d];
73 : 1271872 : const bool not_left_boundary = (global_index != 0);
74 : 1271872 : const bool right_boundary = (global_index == size - 1);
75 : :
76 : 1271872 : coords[d] -= 1;
77 [ + - ]: 1271872 : auto&& left = apply(u_m, coords);
78 : :
79 : 1271872 : coords[d] += 2;
80 [ + - ]: 1271872 : auto&& right = apply(u_m, coords);
81 : :
82 : : // not_left_boundary and right_boundary are boolean values
83 : : // Because of periodic boundary conditions we need to add this boolean mask to
84 : : // obtain the lower triangular part of the Laplace Operator
85 : 1271872 : res += hvector_m[d] * (not_left_boundary * left + right_boundary * right);
86 : : }
87 : 238656 : return res;
88 : : }
89 : :
90 : : private:
91 : : using Mesh_t = typename E::Mesh_t;
92 : : using Layout_t = typename E::Layout_t;
93 : : using vector_type = typename Mesh_t::vector_type;
94 : : using domain_type = typename Layout_t::NDIndex_t;
95 : : const E u_m;
96 : : const vector_type hvector_m;
97 : : const unsigned nghosts_m;
98 : : const domain_type ldom_m;
99 : : const domain_type domain_m;
100 : : };
101 : :
102 : : template <typename E>
103 : : struct meta_upper_laplace
104 : : : public Expression<meta_upper_laplace<E>,
105 : : sizeof(E) + sizeof(typename E::Mesh_t::vector_type)
106 : : + 2 * sizeof(typename E::Layout_t::NDIndex_t)
107 : : + sizeof(unsigned)> {
108 : : constexpr static unsigned dim = E::dim;
109 : :
110 : : KOKKOS_FUNCTION
111 : 12 : meta_upper_laplace(const E& u, const typename E::Mesh_t::vector_type& hvector,
112 : : unsigned nghosts, const typename E::Layout_t::NDIndex_t& ldom,
113 : : const typename E::Layout_t::NDIndex_t& domain)
114 : 12 : : u_m(u)
115 : 12 : , hvector_m(hvector)
116 : 12 : , nghosts_m(nghosts)
117 : 12 : , ldom_m(ldom)
118 : 12 : , domain_m(domain) {}
119 : :
120 : : /*
121 : : * n-dimensional upper triangular Laplacian
122 : : */
123 : : template <typename... Idx>
124 : 238656 : KOKKOS_INLINE_FUNCTION auto operator()(const Idx... args) const {
125 : : using index_type = std::tuple_element_t<0, std::tuple<Idx...>>;
126 : : using T = typename E::Mesh_t::value_type;
127 : 238656 : T res = 0;
128 : :
129 [ + + ]: 1510528 : for (unsigned d = 0; d < dim; d++) {
130 : 1271872 : index_type coords[dim] = {args...};
131 : 1271872 : const int global_index = coords[d] + ldom_m[d].first() - nghosts_m;
132 [ + - ]: 1271872 : const int size = domain_m.length()[d];
133 : 1271872 : const bool left_boundary = (global_index == 0);
134 : 1271872 : const bool not_right_boundary = (global_index != size - 1);
135 : :
136 : 1271872 : coords[d] -= 1;
137 [ + - ]: 1271872 : auto&& left = apply(u_m, coords);
138 : :
139 : 1271872 : coords[d] += 2;
140 [ + - ]: 1271872 : auto&& right = apply(u_m, coords);
141 : :
142 : : // left_boundary and not_right_boundary are boolean values
143 : : // Because of periodic boundary conditions we need to add this boolean mask to
144 : : // obtain the upper triangular part of the Laplace Operator
145 : 1271872 : res += hvector_m[d] * (left_boundary * left + not_right_boundary * right);
146 : : }
147 : 238656 : return res;
148 : : }
149 : :
150 : : private:
151 : : using Mesh_t = typename E::Mesh_t;
152 : : using Layout_t = typename E::Layout_t;
153 : : using vector_type = typename Mesh_t::vector_type;
154 : : using domain_type = typename Layout_t::NDIndex_t;
155 : : const E u_m;
156 : : const vector_type hvector_m;
157 : : const unsigned nghosts_m;
158 : : const domain_type ldom_m;
159 : : const domain_type domain_m;
160 : : };
161 : :
162 : : template <typename E>
163 : : struct meta_upper_and_lower_laplace
164 : : : public Expression<meta_upper_and_lower_laplace<E>,
165 : : sizeof(E) + sizeof(typename E::Mesh_t::vector_type)> {
166 : : constexpr static unsigned dim = E::dim;
167 : :
168 : : KOKKOS_FUNCTION
169 : 12 : meta_upper_and_lower_laplace(const E& u, const typename E::Mesh_t::vector_type& hvector)
170 : 12 : : u_m(u)
171 : 12 : , hvector_m(hvector) {}
172 : :
173 : : /*
174 : : * n-dimensional upper+lower triangular Laplacian
175 : : */
176 : : template <typename... Idx>
177 : 238656 : KOKKOS_INLINE_FUNCTION auto operator()(const Idx... args) const {
178 : : using index_type = std::tuple_element_t<0, std::tuple<Idx...>>;
179 : : using T = typename E::Mesh_t::value_type;
180 : 238656 : T res = 0;
181 [ + + ]: 1510528 : for (unsigned d = 0; d < dim; d++) {
182 : 1271872 : index_type coords[dim] = {args...};
183 : 1271872 : coords[d] -= 1;
184 [ + - ]: 1271872 : auto&& left = apply(u_m, coords);
185 : 1271872 : coords[d] += 2;
186 [ + - ]: 1271872 : auto&& right = apply(u_m, coords);
187 : 1271872 : res += hvector_m[d] * (left + right);
188 : : }
189 : 238656 : return res;
190 : : }
191 : :
192 : : private:
193 : : using vector_type = typename E::Mesh_t::vector_type;
194 : : const E u_m;
195 : : const vector_type hvector_m;
196 : : };
197 : : } // namespace detail
198 : :
199 : : /*!
200 : : * User interface of poisson
201 : : * @param u field
202 : : */
203 : : template <typename Field>
204 : : detail::meta_poisson<Field> poisson(Field& u) {
205 : : constexpr unsigned Dim = Field::dim;
206 : :
207 : : u.fillHalo();
208 : : BConds<Field, Dim>& bcField = u.getFieldBC();
209 : : bcField.apply(u);
210 : :
211 : : return detail::meta_poisson<Field>(u);
212 : : }
213 : :
214 : : /*!
215 : : * User interface of lower triangular Laplacian
216 : : * @param u field
217 : : */
218 : : template <typename Field>
219 : 12 : detail::meta_lower_laplace<Field> lower_laplace(Field& u) {
220 : 12 : constexpr unsigned Dim = Field::dim;
221 : :
222 : 12 : u.fillHalo();
223 : 12 : BConds<Field, Dim>& bcField = u.getFieldBC();
224 : 12 : bcField.apply(u);
225 : :
226 : 12 : return lower_laplace_no_comm(u);
227 : : }
228 : :
229 : : /*!
230 : : * User interface of lower triangular Laplacian without exchange of halo cells
231 : : * @param u field
232 : : */
233 : : template <typename Field>
234 : 12 : detail::meta_lower_laplace<Field> lower_laplace_no_comm(Field& u) {
235 : 12 : constexpr unsigned Dim = Field::dim;
236 : :
237 : : using mesh_type = typename Field::Mesh_t;
238 : 12 : mesh_type& mesh = u.get_mesh();
239 [ + - ]: 12 : typename mesh_type::vector_type hvector(0);
240 [ + + ]: 54 : for (unsigned d = 0; d < Dim; d++) {
241 [ + - ]: 42 : hvector[d] = 1.0 / Kokkos::pow(mesh.getMeshSpacing(d), 2);
242 : : }
243 [ + - ]: 12 : const auto& layout = u.getLayout();
244 : 12 : unsigned nghosts = u.getNghost();
245 [ + - ]: 12 : const auto& ldom = layout.getLocalNDIndex();
246 : 12 : const auto& domain = layout.getDomain();
247 [ + - ]: 24 : return detail::meta_lower_laplace<Field>(u, hvector, nghosts, ldom, domain);
248 : 12 : }
249 : :
250 : : /*!
251 : : * User interface of upper triangular Laplacian
252 : : * @param u field
253 : : */
254 : : template <typename Field>
255 : 12 : detail::meta_upper_laplace<Field> upper_laplace(Field& u) {
256 : 12 : constexpr unsigned Dim = Field::dim;
257 : :
258 : 12 : u.fillHalo();
259 : 12 : BConds<Field, Dim>& bcField = u.getFieldBC();
260 : 12 : bcField.apply(u);
261 : :
262 : 12 : return upper_laplace_no_comm(u);
263 : : }
264 : :
265 : : /*!
266 : : * User interface of upper triangular Laplacian without exchange of halo cells
267 : : * @param u field
268 : : */
269 : : template <typename Field>
270 : 12 : detail::meta_upper_laplace<Field> upper_laplace_no_comm(Field& u) {
271 : 12 : constexpr unsigned Dim = Field::dim;
272 : :
273 : : using mesh_type = typename Field::Mesh_t;
274 : 12 : mesh_type& mesh = u.get_mesh();
275 [ + - ]: 12 : typename mesh_type::vector_type hvector(0);
276 [ + + ]: 54 : for (unsigned d = 0; d < Dim; d++) {
277 [ + - ]: 42 : hvector[d] = 1.0 / Kokkos::pow(mesh.getMeshSpacing(d), 2);
278 : : }
279 [ + - ]: 12 : const auto& layout = u.getLayout();
280 : 12 : unsigned nghosts = u.getNghost();
281 [ + - ]: 12 : const auto& ldom = layout.getLocalNDIndex();
282 : 12 : const auto& domain = layout.getDomain();
283 [ + - ]: 24 : return detail::meta_upper_laplace<Field>(u, hvector, nghosts, ldom, domain);
284 : 12 : }
285 : :
286 : : /*!
287 : : * User interface of upper+lower triangular Laplacian
288 : : * @param u field
289 : : */
290 : : template <typename Field>
291 : 12 : detail::meta_upper_and_lower_laplace<Field> upper_and_lower_laplace(Field& u) {
292 : 12 : constexpr unsigned Dim = Field::dim;
293 : :
294 : 12 : u.fillHalo();
295 : 12 : BConds<Field, Dim>& bcField = u.getFieldBC();
296 : 12 : bcField.apply(u);
297 : :
298 : 12 : return upper_and_lower_laplace_no_comm(u);
299 : : }
300 : :
301 : : /*!
302 : : * User interface of upper+lower triangular Laplacian without exchange of halo cells
303 : : * @param u field
304 : : */
305 : : template <typename Field>
306 : 12 : detail::meta_upper_and_lower_laplace<Field> upper_and_lower_laplace_no_comm(Field& u) {
307 : 12 : constexpr unsigned Dim = Field::dim;
308 : :
309 : : using mesh_type = typename Field::Mesh_t;
310 : 12 : mesh_type& mesh = u.get_mesh();
311 [ + - ]: 12 : typename mesh_type::vector_type hvector(0);
312 [ + + ]: 54 : for (unsigned d = 0; d < Dim; d++) {
313 [ + - ]: 42 : hvector[d] = 1.0 / Kokkos::pow(mesh.getMeshSpacing(d), 2);
314 : : }
315 [ + - ]: 24 : return detail::meta_upper_and_lower_laplace<Field>(u, hvector);
316 : 12 : }
317 : :
318 : : /*!
319 : : * Returns the factor by which to multiply the u field to get
320 : : * the inverse of the diagonal of the Laplacian.
321 : : *
322 : : * @param u field
323 : : */
324 : :
325 : : template <typename Field>
326 : 0 : double negative_inverse_diagonal_laplace(Field& u) {
327 : 0 : constexpr unsigned Dim = Field::dim;
328 : : using mesh_type = typename Field::Mesh_t;
329 : 0 : mesh_type& mesh = u.get_mesh();
330 : 0 : double sum = 0.0;
331 : 0 : double factor = 1.0;
332 : 0 : typename mesh_type::vector_type hvector(0);
333 [ # # ]: 0 : for (unsigned d = 0; d < Dim; ++d) {
334 [ # # ]: 0 : hvector[d] = Kokkos::pow(mesh.getMeshSpacing(d), 2);
335 [ # # ]: 0 : sum += hvector[d] * Kokkos::pow(mesh.getMeshSpacing((d + 1) % Dim), 2);
336 : 0 : factor *= hvector[d];
337 : : }
338 : :
339 : 0 : return 0.5 * (factor / sum);
340 : 0 : }
341 : :
342 : : template <typename Field>
343 : 0 : double diagonal_laplace(Field& u) {
344 : 0 : constexpr unsigned Dim = Field::dim;
345 : : using mesh_type = typename Field::Mesh_t;
346 : 0 : mesh_type& mesh = u.get_mesh();
347 : 0 : double sum = 0.0;
348 [ # # ]: 0 : for (unsigned d = 0; d < Dim; ++d) {
349 : 0 : sum += 1 / (Kokkos::pow(mesh.getMeshSpacing(d), 2));
350 : : }
351 : :
352 : : // u = - 2.0 * sum * u;
353 : 0 : return -2.0 * sum;
354 : : }
355 : :
356 : : template <typename Field>
357 : : void mult(Field& u, const double c) {
358 : : using view_type = Field::view_type;
359 : :
360 : : view_type view = u.getView();
361 : :
362 : : Kokkos::parallel_for(
363 : : "Field_mult_const", u.getFieldRangePolicy(),
364 : : KOKKOS_LAMBDA(int i, int j, int k) { view(i, j, k) *= c; });
365 : : return;
366 : : }
367 : : } // namespace ippl
368 : : #endif // IPPL_LAPLACE_HELPERS_H
|
