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