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 24 : 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 24 : : u_m(u)
56 24 : , hvector_m(hvector)
57 24 : , nghosts_m(nghosts)
58 24 : , ldom_m(ldom)
59 24 : , 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 24 : 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 24 : : u_m(u)
117 24 : , hvector_m(hvector)
118 24 : , nghosts_m(nghosts)
119 24 : , ldom_m(ldom)
120 24 : , 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 24 : meta_upper_and_lower_laplace(const E& u, const typename E::Mesh_t::vector_type& hvector)
173 24 : : u_m(u)
174 24 : , 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 24 : detail::meta_lower_laplace<Field> lower_laplace(Field& u) {
223 24 : constexpr unsigned Dim = Field::dim;
224 :
225 24 : u.fillHalo();
226 24 : BConds<Field, Dim>& bcField = u.getFieldBC();
227 24 : bcField.apply(u);
228 :
229 24 : 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 24 : detail::meta_lower_laplace<Field> lower_laplace_no_comm(Field& u) {
238 24 : constexpr unsigned Dim = Field::dim;
239 :
240 : using mesh_type = typename Field::Mesh_t;
241 24 : mesh_type& mesh = u.get_mesh();
242 24 : typename mesh_type::vector_type hvector(0);
243 108 : for (unsigned d = 0; d < Dim; d++) {
244 84 : hvector[d] = 1.0 / Kokkos::pow(mesh.getMeshSpacing(d), 2);
245 : }
246 24 : const auto& layout = u.getLayout();
247 24 : unsigned nghosts = u.getNghost();
248 24 : const auto& ldom = layout.getLocalNDIndex();
249 24 : const auto& domain = layout.getDomain();
250 48 : return detail::meta_lower_laplace<Field>(u, hvector, nghosts, ldom, domain);
251 24 : }
252 :
253 : /*!
254 : * User interface of upper triangular Laplacian
255 : * @param u field
256 : */
257 : template <typename Field>
258 24 : detail::meta_upper_laplace<Field> upper_laplace(Field& u) {
259 24 : constexpr unsigned Dim = Field::dim;
260 :
261 24 : u.fillHalo();
262 24 : BConds<Field, Dim>& bcField = u.getFieldBC();
263 24 : bcField.apply(u);
264 :
265 24 : 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 24 : detail::meta_upper_laplace<Field> upper_laplace_no_comm(Field& u) {
274 24 : constexpr unsigned Dim = Field::dim;
275 :
276 : using mesh_type = typename Field::Mesh_t;
277 24 : mesh_type& mesh = u.get_mesh();
278 24 : typename mesh_type::vector_type hvector(0);
279 108 : for (unsigned d = 0; d < Dim; d++) {
280 84 : hvector[d] = 1.0 / Kokkos::pow(mesh.getMeshSpacing(d), 2);
281 : }
282 24 : const auto& layout = u.getLayout();
283 24 : unsigned nghosts = u.getNghost();
284 24 : const auto& ldom = layout.getLocalNDIndex();
285 24 : const auto& domain = layout.getDomain();
286 48 : return detail::meta_upper_laplace<Field>(u, hvector, nghosts, ldom, domain);
287 24 : }
288 :
289 : /*!
290 : * User interface of upper+lower triangular Laplacian
291 : * @param u field
292 : */
293 : template <typename Field>
294 24 : detail::meta_upper_and_lower_laplace<Field> upper_and_lower_laplace(Field& u) {
295 24 : constexpr unsigned Dim = Field::dim;
296 :
297 24 : u.fillHalo();
298 24 : BConds<Field, Dim>& bcField = u.getFieldBC();
299 24 : bcField.apply(u);
300 :
301 24 : 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 24 : detail::meta_upper_and_lower_laplace<Field> upper_and_lower_laplace_no_comm(Field& u) {
310 24 : constexpr unsigned Dim = Field::dim;
311 :
312 : using mesh_type = typename Field::Mesh_t;
313 24 : mesh_type& mesh = u.get_mesh();
314 24 : typename mesh_type::vector_type hvector(0);
315 108 : for (unsigned d = 0; d < Dim; d++) {
316 84 : hvector[d] = 1.0 / Kokkos::pow(mesh.getMeshSpacing(d), 2);
317 : }
318 48 : return detail::meta_upper_and_lower_laplace<Field>(u, hvector);
319 24 : }
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
|
