Line data Source code
1 : //
2 : // Class FFTOpenPoissonSolver
3 : // FFT-based Poisson Solver for open boundaries.
4 : // Solves laplace(phi) = -rho, and E = -grad(phi).
5 : //
6 : //
7 :
8 : // Communication specific functions (pack and unpack).
9 : template <typename Tb, typename Tf>
10 0 : void pack(const ippl::NDIndex<3> intersect, Kokkos::View<Tf***>& view,
11 : ippl::detail::FieldBufferData<Tb>& fd, int nghost, const ippl::NDIndex<3> ldom,
12 : ippl::mpi::Communicator::size_type& nsends) {
13 0 : Kokkos::View<Tb*>& buffer = fd.buffer;
14 :
15 0 : size_t size = intersect.size();
16 0 : nsends = size;
17 0 : if (buffer.size() < size) {
18 0 : const int overalloc = ippl::Comm->getDefaultOverallocation();
19 0 : Kokkos::realloc(buffer, size * overalloc);
20 : }
21 :
22 0 : const int first0 = intersect[0].first() + nghost - ldom[0].first();
23 0 : const int first1 = intersect[1].first() + nghost - ldom[1].first();
24 0 : const int first2 = intersect[2].first() + nghost - ldom[2].first();
25 :
26 0 : const int last0 = intersect[0].last() + nghost - ldom[0].first() + 1;
27 0 : const int last1 = intersect[1].last() + nghost - ldom[1].first() + 1;
28 0 : const int last2 = intersect[2].last() + nghost - ldom[2].first() + 1;
29 :
30 : using mdrange_type = Kokkos::MDRangePolicy<Kokkos::Rank<3>>;
31 : // This type casting to long int is necessary as otherwise Kokkos complains for
32 : // intel compilers
33 0 : Kokkos::parallel_for(
34 : "pack()",
35 0 : mdrange_type({first0, first1, first2}, {(long int)last0, (long int)last1, (long int)last2}),
36 0 : KOKKOS_LAMBDA(const size_t i, const size_t j, const size_t k) {
37 0 : const int ig = i - first0;
38 0 : const int jg = j - first1;
39 0 : const int kg = k - first2;
40 :
41 0 : int l = ig + jg * intersect[0].length()
42 0 : + kg * intersect[1].length() * intersect[0].length();
43 :
44 0 : Kokkos::complex<Tb> val = view(i, j, k);
45 :
46 0 : buffer(l) = Kokkos::real(val);
47 : });
48 0 : Kokkos::fence();
49 0 : }
50 :
51 : template <int tensorRank, typename Tb, typename Tf>
52 0 : void unpack_impl(const ippl::NDIndex<3> intersect, const Kokkos::View<Tf***>& view,
53 : ippl::detail::FieldBufferData<Tb>& fd, int nghost, const ippl::NDIndex<3> ldom,
54 : size_t dim1 = 0, size_t dim2 = 0, bool x = false, bool y = false, bool z = false) {
55 0 : Kokkos::View<Tb*>& buffer = fd.buffer;
56 :
57 0 : const int first0 = intersect[0].first() + nghost - ldom[0].first();
58 0 : const int first1 = intersect[1].first() + nghost - ldom[1].first();
59 0 : const int first2 = intersect[2].first() + nghost - ldom[2].first();
60 :
61 0 : const int last0 = intersect[0].last() + nghost - ldom[0].first() + 1;
62 0 : const int last1 = intersect[1].last() + nghost - ldom[1].first() + 1;
63 0 : const int last2 = intersect[2].last() + nghost - ldom[2].first() + 1;
64 :
65 : using mdrange_type = Kokkos::MDRangePolicy<Kokkos::Rank<3>>;
66 0 : Kokkos::parallel_for(
67 0 : "pack()", mdrange_type({first0, first1, first2}, {last0, last1, last2}),
68 0 : KOKKOS_LAMBDA(const size_t i, const size_t j, const size_t k) {
69 0 : int ig = i - first0;
70 0 : int jg = j - first1;
71 0 : int kg = k - first2;
72 :
73 0 : ig = x * (intersect[0].length() - 2 * ig - 1) + ig;
74 0 : jg = y * (intersect[1].length() - 2 * jg - 1) + jg;
75 0 : kg = z * (intersect[2].length() - 2 * kg - 1) + kg;
76 :
77 0 : int l = ig + jg * intersect[0].length()
78 0 : + kg * intersect[1].length() * intersect[0].length();
79 :
80 0 : ippl::detail::ViewAccess<tensorRank, decltype(view)>::get(view, dim1, dim2, i, j, k) =
81 0 : buffer(l);
82 : });
83 0 : Kokkos::fence();
84 0 : }
85 :
86 : template <typename Tb, typename Tf>
87 0 : void unpack(const ippl::NDIndex<3> intersect, const Kokkos::View<Tf***>& view,
88 : ippl::detail::FieldBufferData<Tb>& fd, int nghost, const ippl::NDIndex<3> ldom,
89 : bool x = false, bool y = false, bool z = false) {
90 0 : unpack_impl<0, Tb, Tf>(intersect, view, fd, nghost, ldom, 0, 0, x, y, z);
91 0 : }
92 :
93 : template <typename Tb, typename Tf>
94 0 : void unpack(const ippl::NDIndex<3> intersect, const Kokkos::View<ippl::Vector<Tf, 3>***>& view,
95 : size_t dim1, ippl::detail::FieldBufferData<Tb>& fd, int nghost,
96 : const ippl::NDIndex<3> ldom) {
97 0 : unpack_impl<1, Tb, ippl::Vector<Tf, 3>>(intersect, view, fd, nghost, ldom, dim1);
98 0 : }
99 :
100 : template <typename Tb, typename Tf>
101 0 : void unpack(const ippl::NDIndex<3> intersect,
102 : const Kokkos::View<ippl::Vector<ippl::Vector<Tf, 3>, 3>***>& view,
103 : ippl::detail::FieldBufferData<Tb>& fd, int nghost, const ippl::NDIndex<3> ldom,
104 : size_t dim1, size_t dim2) {
105 0 : unpack_impl<2, Tb, ippl::Vector<ippl::Vector<Tf, 3>, 3>>(intersect, view, fd, nghost, ldom,
106 : dim1, dim2);
107 0 : }
108 :
109 : template <typename Tb, typename Tf, unsigned Dim>
110 0 : void solver_send(int TAG, int id, int i, const ippl::NDIndex<Dim> intersection,
111 : const ippl::NDIndex<Dim> ldom, int nghost, Kokkos::View<Tf***>& view,
112 : ippl::detail::FieldBufferData<Tb>& fd, std::vector<MPI_Request>& requests) {
113 : using memory_space = typename Kokkos::View<Tf***>::memory_space;
114 :
115 0 : requests.resize(requests.size() + 1);
116 :
117 : ippl::mpi::Communicator::size_type nsends;
118 0 : pack(intersection, view, fd, nghost, ldom, nsends);
119 :
120 : // Buffer message indicates the domain intersection (x, y, z, xy, yz, xz, xyz).
121 0 : ippl::mpi::Communicator::buffer_type<memory_space> buf =
122 : ippl::Comm->getBuffer<memory_space, Tf>(nsends);
123 :
124 0 : int tag = TAG + id;
125 :
126 0 : ippl::Comm->isend(i, tag, fd, *buf, requests.back(), nsends);
127 0 : buf->resetWritePos();
128 0 : }
129 :
130 : template <typename Tb, typename Tf, unsigned Dim>
131 0 : void solver_recv(int TAG, int id, int i, const ippl::NDIndex<Dim> intersection,
132 : const ippl::NDIndex<Dim> ldom, int nghost, Kokkos::View<Tf***>& view,
133 : ippl::detail::FieldBufferData<Tb>& fd, bool x = false, bool y = false,
134 : bool z = false) {
135 : using memory_space = typename Kokkos::View<Tf***>::memory_space;
136 :
137 : ippl::mpi::Communicator::size_type nrecvs;
138 0 : nrecvs = intersection.size();
139 :
140 : // Buffer message indicates the domain intersection (x, y, z, xy, yz, xz, xyz).
141 0 : ippl::mpi::Communicator::buffer_type<memory_space> buf =
142 : ippl::Comm->getBuffer<memory_space, Tf>(nrecvs);
143 :
144 0 : int tag = TAG + id;
145 :
146 0 : ippl::Comm->recv(i, tag, fd, *buf, nrecvs * sizeof(Tf), nrecvs);
147 0 : buf->resetReadPos();
148 :
149 0 : unpack(intersection, view, fd, nghost, ldom, x, y, z);
150 0 : }
151 :
152 : namespace ippl {
153 :
154 : /////////////////////////////////////////////////////////////////////////
155 : // constructor and destructor
156 : template <typename FieldLHS, typename FieldRHS>
157 : FFTOpenPoissonSolver<FieldLHS, FieldRHS>::FFTOpenPoissonSolver()
158 : : Base()
159 : , mesh_mp(nullptr)
160 : , layout_mp(nullptr)
161 : , mesh2_m(nullptr)
162 : , layout2_m(nullptr)
163 : , meshComplex_m(nullptr)
164 : , layoutComplex_m(nullptr)
165 : , mesh4_m(nullptr)
166 : , layout4_m(nullptr)
167 : , mesh2n1_m(nullptr)
168 : , layout2n1_m(nullptr)
169 : , isGradFD_m(false) {
170 : setDefaultParameters();
171 : }
172 :
173 : template <typename FieldLHS, typename FieldRHS>
174 0 : FFTOpenPoissonSolver<FieldLHS, FieldRHS>::FFTOpenPoissonSolver(rhs_type& rhs,
175 : ParameterList& params)
176 0 : : mesh_mp(nullptr)
177 0 : , layout_mp(nullptr)
178 0 : , mesh2_m(nullptr)
179 0 : , layout2_m(nullptr)
180 0 : , meshComplex_m(nullptr)
181 0 : , layoutComplex_m(nullptr)
182 0 : , mesh4_m(nullptr)
183 0 : , layout4_m(nullptr)
184 0 : , mesh2n1_m(nullptr)
185 0 : , layout2n1_m(nullptr)
186 0 : , isGradFD_m(false) {
187 : using T = typename FieldLHS::value_type::value_type;
188 : static_assert(std::is_floating_point<T>::value, "Not a floating point type");
189 :
190 0 : setDefaultParameters();
191 0 : this->params_m.merge(params);
192 0 : this->params_m.update("output_type", Base::SOL);
193 :
194 0 : this->setRhs(rhs);
195 0 : }
196 :
197 : template <typename FieldLHS, typename FieldRHS>
198 0 : FFTOpenPoissonSolver<FieldLHS, FieldRHS>::FFTOpenPoissonSolver(lhs_type& lhs, rhs_type& rhs,
199 : ParameterList& params)
200 0 : : mesh_mp(nullptr)
201 0 : , layout_mp(nullptr)
202 0 : , mesh2_m(nullptr)
203 0 : , layout2_m(nullptr)
204 0 : , meshComplex_m(nullptr)
205 0 : , layoutComplex_m(nullptr)
206 0 : , mesh4_m(nullptr)
207 0 : , layout4_m(nullptr)
208 0 : , mesh2n1_m(nullptr)
209 0 : , layout2n1_m(nullptr)
210 0 : , isGradFD_m(false) {
211 : using T = typename FieldLHS::value_type::value_type;
212 : static_assert(std::is_floating_point<T>::value, "Not a floating point type");
213 :
214 0 : setDefaultParameters();
215 0 : this->params_m.merge(params);
216 :
217 0 : this->setLhs(lhs);
218 0 : this->setRhs(rhs);
219 0 : }
220 :
221 : /////////////////////////////////////////////////////////////////////////
222 : // override setRhs to call class-specific initialization
223 : template <typename FieldLHS, typename FieldRHS>
224 0 : void FFTOpenPoissonSolver<FieldLHS, FieldRHS>::setRhs(rhs_type& rhs) {
225 0 : Base::setRhs(rhs);
226 :
227 : // start a timer
228 0 : static IpplTimings::TimerRef initialize = IpplTimings::getTimer("Initialize");
229 0 : IpplTimings::startTimer(initialize);
230 :
231 0 : initializeFields();
232 :
233 0 : IpplTimings::stopTimer(initialize);
234 0 : }
235 :
236 : /////////////////////////////////////////////////////////////////////////
237 : // allows user to set gradient of phi = Efield instead of spectral
238 : // calculation of Efield (which uses FFTs)
239 :
240 : template <typename FieldLHS, typename FieldRHS>
241 : void FFTOpenPoissonSolver<FieldLHS, FieldRHS>::setGradFD() {
242 : // get the output type (sol, grad, or sol & grad)
243 : const int out = this->params_m.template get<int>("output_type");
244 :
245 : if (out != Base::SOL_AND_GRAD) {
246 : throw IpplException(
247 : "FFTOpenPoissonSolver::setGradFD()",
248 : "Cannot use gradient for Efield computation unless output type is SOL_AND_GRAD");
249 : } else {
250 : isGradFD_m = true;
251 : }
252 : }
253 :
254 : /////////////////////////////////////////////////////////////////////////
255 : // initializeFields method, called in constructor
256 :
257 : template <typename FieldLHS, typename FieldRHS>
258 0 : void FFTOpenPoissonSolver<FieldLHS, FieldRHS>::initializeFields() {
259 : // get algorithm and hessian flag from parameter list
260 0 : const int alg = this->params_m.template get<int>("algorithm");
261 0 : const bool hessian = this->params_m.template get<bool>("hessian");
262 :
263 : // first check if valid algorithm choice
264 0 : if ((alg != Algorithm::VICO) && (alg != Algorithm::HOCKNEY)
265 0 : && (alg != Algorithm::BIHARMONIC) && (alg != Algorithm::DCT_VICO)) {
266 0 : throw IpplException("FFTOpenPoissonSolver::initializeFields()",
267 : "Currently only HOCKNEY, VICO, DCT_VICO, and BIHARMONIC are "
268 : "supported for open BCs");
269 : }
270 :
271 : // get layout and mesh
272 0 : layout_mp = &(this->rhs_mp->getLayout());
273 0 : mesh_mp = &(this->rhs_mp->get_mesh());
274 0 : mpi::Communicator comm = layout_mp->comm;
275 :
276 : // get mesh spacing and origin
277 0 : hr_m = mesh_mp->getMeshSpacing();
278 0 : vector_type origin = mesh_mp->getOrigin();
279 :
280 : // create domain for the real fields
281 0 : domain_m = layout_mp->getDomain();
282 :
283 : // get the mesh spacings and sizes for each dimension
284 0 : for (unsigned int i = 0; i < Dim; ++i) {
285 0 : nr_m[i] = domain_m[i].length();
286 :
287 : // create the doubled domain for the FFT procedure
288 0 : domain2_m[i] = Index(2 * nr_m[i]);
289 : }
290 :
291 : // define decomposition (parallel / serial)
292 0 : std::array<bool, Dim> isParallel = layout_mp->isParallel();
293 :
294 : // create double sized mesh and layout objects using the previously defined domain2_m
295 : using mesh_type = typename lhs_type::Mesh_t;
296 0 : mesh2_m = std::unique_ptr<mesh_type>(new mesh_type(domain2_m, hr_m, origin));
297 0 : layout2_m = std::unique_ptr<FieldLayout_t>(new FieldLayout_t(comm, domain2_m, isParallel));
298 :
299 : // create the domain for the transformed (complex) fields
300 : // since we use HeFFTe for the transforms it doesn't require permuting to the right
301 : // one of the dimensions has only (n/2 +1) as our original fields are fully real
302 : // the dimension is given by the user via r2c_direction
303 0 : unsigned int RCDirection = this->params_m.template get<int>("r2c_direction");
304 0 : for (unsigned int i = 0; i < Dim; ++i) {
305 0 : if (i == RCDirection) {
306 0 : domainComplex_m[RCDirection] = Index(nr_m[RCDirection] + 1);
307 : } else {
308 0 : domainComplex_m[i] = Index(2 * nr_m[i]);
309 : }
310 : }
311 :
312 : // create mesh and layout for the real to complex FFT transformed fields
313 0 : meshComplex_m = std::unique_ptr<mesh_type>(new mesh_type(domainComplex_m, hr_m, origin));
314 0 : layoutComplex_m =
315 0 : std::unique_ptr<FieldLayout_t>(new FieldLayout_t(comm, domainComplex_m, isParallel));
316 :
317 : // initialize fields
318 0 : storage_field.initialize(*mesh2_m, *layout2_m);
319 0 : rho2tr_m.initialize(*meshComplex_m, *layoutComplex_m);
320 0 : grntr_m.initialize(*meshComplex_m, *layoutComplex_m);
321 :
322 0 : int out = this->params_m.template get<int>("output_type");
323 0 : if (((out == Base::GRAD || out == Base::SOL_AND_GRAD) && !isGradFD_m) || hessian) {
324 0 : temp_m.initialize(*meshComplex_m, *layoutComplex_m);
325 : }
326 :
327 0 : if (hessian) {
328 0 : hess_m.initialize(*mesh_mp, *layout_mp);
329 : }
330 :
331 : // create the FFT object
332 0 : fft_m = std::make_unique<FFT_t>(*layout2_m, *layoutComplex_m, this->params_m);
333 :
334 : // if Vico, also need to create mesh and layout for 4N Fourier domain
335 : // on this domain, the truncated Green's function is defined
336 : // also need to create the 4N complex grid, on which precomputation step done
337 0 : if (alg == Algorithm::VICO || alg == Algorithm::BIHARMONIC) {
338 : // start a timer
339 0 : static IpplTimings::TimerRef initialize_vico =
340 0 : IpplTimings::getTimer("Initialize: extra Vico");
341 0 : IpplTimings::startTimer(initialize_vico);
342 :
343 0 : for (unsigned int i = 0; i < Dim; ++i) {
344 0 : domain4_m[i] = Index(4 * nr_m[i]);
345 : }
346 :
347 : // 4N grid
348 : using mesh_type = typename lhs_type::Mesh_t;
349 0 : mesh4_m = std::unique_ptr<mesh_type>(new mesh_type(domain4_m, hr_m, origin));
350 0 : layout4_m =
351 0 : std::unique_ptr<FieldLayout_t>(new FieldLayout_t(comm, domain4_m, isParallel));
352 :
353 : // initialize fields
354 0 : grnL_m.initialize(*mesh4_m, *layout4_m);
355 :
356 : // create a Complex-to-Complex FFT object to transform for layout4
357 0 : fft4n_m = std::make_unique<FFT<CCTransform, CxField_gt>>(*layout4_m, this->params_m);
358 :
359 0 : IpplTimings::stopTimer(initialize_vico);
360 : }
361 :
362 : // if DCT_VICO, need 2N+1 mesh, layout, domain, and green's function field for
363 : // precomputation
364 0 : if (alg == Algorithm::DCT_VICO) {
365 : // start a timer
366 0 : static IpplTimings::TimerRef initialize_vico =
367 0 : IpplTimings::getTimer("Initialize: extra Vico");
368 0 : IpplTimings::startTimer(initialize_vico);
369 :
370 : // 2N+1 domain for DCT
371 0 : for (unsigned int i = 0; i < Dim; ++i) {
372 0 : domain2n1_m[i] = Index(2 * nr_m[i] + 1);
373 : }
374 :
375 : // 2N+1 grid
376 0 : mesh2n1_m = std::unique_ptr<mesh_type>(new mesh_type(domain2n1_m, hr_m, origin));
377 0 : layout2n1_m =
378 0 : std::unique_ptr<FieldLayout_t>(new FieldLayout_t(comm, domain2n1_m, isParallel));
379 :
380 : // initialize fields
381 0 : grn2n1_m.initialize(*mesh2n1_m, *layout2n1_m); // 2N+1 grnL
382 :
383 : // create real to real FFT object for 2N+1 grid
384 0 : fft2n1_m = std::make_unique<FFT<Cos1Transform, Field_t>>(*layout2n1_m, this->params_m);
385 :
386 0 : IpplTimings::stopTimer(initialize_vico);
387 : }
388 :
389 : // these are fields that are used for calculating the Green's function for Hockney
390 0 : if (alg == Algorithm::HOCKNEY) {
391 : // start a timer
392 0 : static IpplTimings::TimerRef initialize_hockney =
393 0 : IpplTimings::getTimer("Initialize: extra Hockney");
394 0 : IpplTimings::startTimer(initialize_hockney);
395 :
396 0 : for (unsigned int d = 0; d < Dim; ++d) {
397 0 : grnIField_m[d].initialize(*mesh2_m, *layout2_m);
398 :
399 : // get number of ghost points and the Kokkos view to iterate over field
400 0 : auto view = grnIField_m[d].getView();
401 0 : const int nghost = grnIField_m[d].getNghost();
402 0 : const auto& ldom = layout2_m->getLocalNDIndex();
403 :
404 : // the length of the physical domain
405 0 : const int size = nr_m[d];
406 :
407 : // Kokkos parallel for loop to initialize grnIField[d]
408 0 : switch (d) {
409 0 : case 0:
410 0 : Kokkos::parallel_for(
411 : "Helper index Green field initialization",
412 0 : grnIField_m[d].getFieldRangePolicy(),
413 0 : KOKKOS_LAMBDA(const int i, const int j, const int k) {
414 : // go from local indices to global
415 0 : const int ig = i + ldom[0].first() - nghost;
416 0 : const int jg = j + ldom[1].first() - nghost;
417 0 : const int kg = k + ldom[2].first() - nghost;
418 :
419 : // assign (index)^2 if 0 <= index < N, and (2N-index)^2 elsewhere
420 0 : const bool outsideN = (ig >= size);
421 0 : view(i, j, k) =
422 0 : (2 * size * outsideN - ig) * (2 * size * outsideN - ig);
423 :
424 : // add 1.0 if at (0,0,0) to avoid singularity
425 0 : const bool isOrig = ((ig == 0) && (jg == 0) && (kg == 0));
426 0 : view(i, j, k) += isOrig * 1.0;
427 : });
428 0 : break;
429 0 : case 1:
430 0 : Kokkos::parallel_for(
431 : "Helper index Green field initialization",
432 0 : grnIField_m[d].getFieldRangePolicy(),
433 0 : KOKKOS_LAMBDA(const int i, const int j, const int k) {
434 : // go from local indices to global
435 0 : const int jg = j + ldom[1].first() - nghost;
436 :
437 : // assign (index)^2 if 0 <= index < N, and (2N-index)^2 elsewhere
438 0 : const bool outsideN = (jg >= size);
439 0 : view(i, j, k) =
440 0 : (2 * size * outsideN - jg) * (2 * size * outsideN - jg);
441 : });
442 0 : break;
443 0 : case 2:
444 0 : Kokkos::parallel_for(
445 : "Helper index Green field initialization",
446 0 : grnIField_m[d].getFieldRangePolicy(),
447 0 : KOKKOS_LAMBDA(const int i, const int j, const int k) {
448 : // go from local indices to global
449 0 : const int kg = k + ldom[2].first() - nghost;
450 :
451 : // assign (index)^2 if 0 <= index < N, and (2N-index)^2 elsewhere
452 0 : const bool outsideN = (kg >= size);
453 0 : view(i, j, k) =
454 0 : (2 * size * outsideN - kg) * (2 * size * outsideN - kg);
455 : });
456 0 : break;
457 : }
458 : }
459 0 : IpplTimings::stopTimer(initialize_hockney);
460 : }
461 :
462 0 : static IpplTimings::TimerRef warmup = IpplTimings::getTimer("Warmup");
463 0 : IpplTimings::startTimer(warmup);
464 :
465 : // "empty" transforms to warmup all the FFTs
466 0 : fft_m->warmup(rho2_mr, rho2tr_m);
467 0 : if (alg == Algorithm::VICO || alg == Algorithm::BIHARMONIC) {
468 0 : fft4n_m->warmup(grnL_m);
469 : }
470 0 : if (alg == Algorithm::DCT_VICO) {
471 0 : fft2n1_m->warmup(grn2n1_m);
472 : }
473 :
474 0 : IpplTimings::stopTimer(warmup);
475 :
476 0 : rho2_mr = 0.0;
477 0 : rho2tr_m = 0.0;
478 0 : grnL_m = 0.0;
479 0 : grn2n1_m = 0.0;
480 :
481 : // call greensFunction and we will get the transformed G in the class attribute
482 : // this is done in initialization so that we already have the precomputed fct
483 : // for all timesteps (green's fct will only change if mesh size changes)
484 0 : static IpplTimings::TimerRef ginit = IpplTimings::getTimer("Green Init");
485 0 : IpplTimings::startTimer(ginit);
486 0 : greensFunction();
487 0 : IpplTimings::stopTimer(ginit);
488 0 : };
489 :
490 : /////////////////////////////////////////////////////////////////////////
491 : // compute electric potential by solving Poisson's eq given a field rho and mesh spacings hr
492 : template <typename FieldLHS, typename FieldRHS>
493 0 : void FFTOpenPoissonSolver<FieldLHS, FieldRHS>::solve() {
494 : // start a timer
495 0 : static IpplTimings::TimerRef solve = IpplTimings::getTimer("Solve");
496 0 : IpplTimings::startTimer(solve);
497 :
498 : // get the output type (sol, grad, or sol & grad)
499 0 : const int out = this->params_m.template get<int>("output_type");
500 :
501 : // get the algorithm (hockney, vico, or biharmonic)
502 0 : const int alg = this->params_m.template get<int>("algorithm");
503 :
504 : // get hessian flag (if true, we compute the Hessian)
505 0 : const bool hessian = this->params_m.template get<bool>("hessian");
506 :
507 : // set the mesh & spacing, which may change each timestep
508 0 : mesh_mp = &(this->rhs_mp->get_mesh());
509 :
510 : // check whether the mesh spacing has changed with respect to the old one
511 : // if yes, update and set green flag to true
512 0 : bool green = false;
513 0 : for (unsigned int i = 0; i < Dim; ++i) {
514 0 : if (hr_m[i] != mesh_mp->getMeshSpacing(i)) {
515 0 : hr_m[i] = mesh_mp->getMeshSpacing(i);
516 0 : green = true;
517 : }
518 : }
519 :
520 : // set mesh spacing on the other grids again
521 0 : mesh2_m->setMeshSpacing(hr_m);
522 0 : meshComplex_m->setMeshSpacing(hr_m);
523 :
524 : // field object on the doubled grid; zero-padded
525 0 : rho2_mr = 0.0;
526 :
527 : // start a timer
528 0 : static IpplTimings::TimerRef stod = IpplTimings::getTimer("Solve: Physical to double");
529 0 : IpplTimings::startTimer(stod);
530 :
531 : // store rho (RHS) in the lower left quadrant of the doubled grid
532 : // with or without communication (if only 1 rank)
533 :
534 0 : const int ranks = Comm->size();
535 :
536 0 : auto view2 = rho2_mr.getView();
537 0 : auto view1 = this->rhs_mp->getView();
538 :
539 0 : const int nghost2 = rho2_mr.getNghost();
540 0 : const int nghost1 = this->rhs_mp->getNghost();
541 :
542 0 : const auto& ldom2 = layout2_m->getLocalNDIndex();
543 0 : const auto& ldom1 = layout_mp->getLocalNDIndex();
544 :
545 0 : if (ranks > 1) {
546 : // COMMUNICATION
547 0 : const auto& lDomains2 = layout2_m->getHostLocalDomains();
548 :
549 : // send
550 0 : std::vector<MPI_Request> requests(0);
551 :
552 0 : for (int i = 0; i < ranks; ++i) {
553 0 : if (lDomains2[i].touches(ldom1)) {
554 0 : auto intersection = lDomains2[i].intersect(ldom1);
555 :
556 0 : solver_send(mpi::tag::OPEN_SOLVER, 0, i, intersection, ldom1, nghost1, view1,
557 0 : fd_m, requests);
558 : }
559 : }
560 :
561 : // receive
562 0 : const auto& lDomains1 = layout_mp->getHostLocalDomains();
563 :
564 0 : for (int i = 0; i < ranks; ++i) {
565 0 : if (lDomains1[i].touches(ldom2)) {
566 0 : auto intersection = lDomains1[i].intersect(ldom2);
567 :
568 : mpi::Communicator::size_type nrecvs;
569 0 : nrecvs = intersection.size();
570 :
571 0 : buffer_type buf = Comm->getBuffer<memory_space, Trhs>(nrecvs);
572 :
573 0 : Comm->recv(i, mpi::tag::OPEN_SOLVER, fd_m, *buf, nrecvs * sizeof(Trhs), nrecvs);
574 0 : buf->resetReadPos();
575 :
576 0 : unpack(intersection, view2, fd_m, nghost2, ldom2);
577 0 : }
578 : }
579 :
580 : // wait for all messages to be received
581 0 : if (requests.size() > 0) {
582 0 : MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
583 : }
584 0 : ippl::Comm->freeAllBuffers();
585 :
586 0 : } else {
587 0 : Kokkos::parallel_for(
588 0 : "Write rho on the doubled grid", this->rhs_mp->getFieldRangePolicy(),
589 0 : KOKKOS_LAMBDA(const size_t i, const size_t j, const size_t k) {
590 0 : const size_t ig2 = i + ldom2[0].first() - nghost2;
591 0 : const size_t jg2 = j + ldom2[1].first() - nghost2;
592 0 : const size_t kg2 = k + ldom2[2].first() - nghost2;
593 :
594 0 : const size_t ig1 = i + ldom1[0].first() - nghost1;
595 0 : const size_t jg1 = j + ldom1[1].first() - nghost1;
596 0 : const size_t kg1 = k + ldom1[2].first() - nghost1;
597 :
598 : // write physical rho on [0,N-1] of doubled field
599 0 : const bool isQuadrant1 = ((ig1 == ig2) && (jg1 == jg2) && (kg1 == kg2));
600 0 : view2(i, j, k) = view1(i, j, k) * isQuadrant1;
601 : });
602 : }
603 :
604 0 : IpplTimings::stopTimer(stod);
605 :
606 : // start a timer
607 0 : static IpplTimings::TimerRef fftrho = IpplTimings::getTimer("FFT: Rho");
608 0 : IpplTimings::startTimer(fftrho);
609 :
610 : // forward FFT of the charge density field on doubled grid
611 0 : fft_m->transform(FORWARD, rho2_mr, rho2tr_m);
612 :
613 0 : IpplTimings::stopTimer(fftrho);
614 :
615 : // call greensFunction to recompute if the mesh spacing has changed
616 0 : if (green) {
617 0 : greensFunction();
618 : }
619 :
620 : // multiply FFT(rho2)*FFT(green)
621 : // convolution becomes multiplication in FFT
622 : // minus sign since we are solving laplace(phi) = -rho
623 0 : rho2tr_m = -rho2tr_m * grntr_m;
624 :
625 : // if output_type is SOL or SOL_AND_GRAD, we caculate solution
626 0 : if ((out == Base::SOL) || (out == Base::SOL_AND_GRAD)) {
627 : // start a timer
628 0 : static IpplTimings::TimerRef fftc = IpplTimings::getTimer("FFT: Convolution");
629 0 : IpplTimings::startTimer(fftc);
630 :
631 : // inverse FFT of the product and store the electrostatic potential in rho2_mr
632 0 : fft_m->transform(BACKWARD, rho2_mr, rho2tr_m);
633 :
634 0 : IpplTimings::stopTimer(fftc);
635 : // Hockney: multiply the rho2_mr field by the total number of points to account for
636 : // double counting (rho and green) of normalization factor in forward transform
637 : // also multiply by the mesh spacing^3 (to account for discretization)
638 : // Vico: need to multiply by normalization factor of 1/4N^3,
639 : // since only backward transform was performed on the 4N grid
640 : // DCT_VICO: need to multiply by a factor of (2N)^3 to match the normalization factor in
641 : // the transform.
642 0 : for (unsigned int i = 0; i < Dim; ++i) {
643 0 : switch (alg) {
644 0 : case Algorithm::HOCKNEY:
645 0 : rho2_mr = rho2_mr * 2.0 * nr_m[i] * hr_m[i];
646 0 : break;
647 0 : case Algorithm::VICO:
648 : case Algorithm::BIHARMONIC:
649 0 : rho2_mr = rho2_mr * 2.0 * (1.0 / 4.0);
650 0 : break;
651 0 : case Algorithm::DCT_VICO:
652 0 : rho2_mr = rho2_mr * (1.0 / 4.0);
653 0 : break;
654 0 : default:
655 0 : throw IpplException(
656 : "FFTOpenPoissonSolver::initializeFields()",
657 : "Currently only HOCKNEY, VICO, DCT_VICO, and BIHARMONIC are "
658 : "supported for open BCs");
659 : }
660 : }
661 :
662 : // start a timer
663 0 : static IpplTimings::TimerRef dtos = IpplTimings::getTimer("Solve: Double to physical");
664 0 : IpplTimings::startTimer(dtos);
665 :
666 : // get the physical part only --> physical electrostatic potential is now given in RHS
667 : // need communication if more than one rank
668 :
669 0 : if (ranks > 1) {
670 : // COMMUNICATION
671 :
672 : // send
673 0 : const auto& lDomains1 = layout_mp->getHostLocalDomains();
674 :
675 0 : std::vector<MPI_Request> requests(0);
676 :
677 0 : for (int i = 0; i < ranks; ++i) {
678 0 : if (lDomains1[i].touches(ldom2)) {
679 0 : auto intersection = lDomains1[i].intersect(ldom2);
680 :
681 0 : solver_send(mpi::tag::OPEN_SOLVER, 0, i, intersection, ldom2, nghost2,
682 0 : view2, fd_m, requests);
683 : }
684 : }
685 :
686 : // receive
687 0 : const auto& lDomains2 = layout2_m->getHostLocalDomains();
688 :
689 0 : for (int i = 0; i < ranks; ++i) {
690 0 : if (ldom1.touches(lDomains2[i])) {
691 0 : auto intersection = ldom1.intersect(lDomains2[i]);
692 :
693 : mpi::Communicator::size_type nrecvs;
694 0 : nrecvs = intersection.size();
695 :
696 0 : buffer_type buf = Comm->getBuffer<memory_space, Trhs>(nrecvs);
697 :
698 0 : Comm->recv(i, mpi::tag::OPEN_SOLVER, fd_m, *buf, nrecvs * sizeof(Trhs),
699 : nrecvs);
700 0 : buf->resetReadPos();
701 :
702 0 : unpack(intersection, view1, fd_m, nghost1, ldom1);
703 0 : }
704 : }
705 :
706 : // wait for all messages to be received
707 0 : if (requests.size() > 0) {
708 0 : MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
709 : }
710 0 : ippl::Comm->freeAllBuffers();
711 :
712 0 : } else {
713 0 : Kokkos::parallel_for(
714 : "Write the solution into the LHS on physical grid",
715 0 : this->rhs_mp->getFieldRangePolicy(),
716 0 : KOKKOS_LAMBDA(const int i, const int j, const int k) {
717 0 : const int ig2 = i + ldom2[0].first() - nghost2;
718 0 : const int jg2 = j + ldom2[1].first() - nghost2;
719 0 : const int kg2 = k + ldom2[2].first() - nghost2;
720 :
721 0 : const int ig = i + ldom1[0].first() - nghost1;
722 0 : const int jg = j + ldom1[1].first() - nghost1;
723 0 : const int kg = k + ldom1[2].first() - nghost1;
724 :
725 : // take [0,N-1] as physical solution
726 0 : const bool isQuadrant1 = ((ig == ig2) && (jg == jg2) && (kg == kg2));
727 0 : view1(i, j, k) = view2(i, j, k) * isQuadrant1;
728 : });
729 : }
730 0 : IpplTimings::stopTimer(dtos);
731 : }
732 :
733 : // if we want finite differences Efield = -grad(phi)
734 : // instead of computing it in Fourier domain
735 : // this is only possible if SOL_AND_GRAD is the output type
736 0 : if (isGradFD_m && (out == Base::SOL_AND_GRAD)) {
737 0 : *(this->lhs_mp) = -grad(*this->rhs_mp);
738 : }
739 :
740 : // if output_type is GRAD or SOL_AND_GRAD, we calculate E-field (gradient in Fourier domain)
741 0 : if (((out == Base::GRAD) || (out == Base::SOL_AND_GRAD)) && (!isGradFD_m)) {
742 : // start a timer
743 0 : static IpplTimings::TimerRef efield = IpplTimings::getTimer("Solve: Electric field");
744 0 : IpplTimings::startTimer(efield);
745 :
746 : // get E field view (LHS)
747 0 : auto viewL = this->lhs_mp->getView();
748 0 : const int nghostL = this->lhs_mp->getNghost();
749 :
750 : // get rho2tr_m view (as we want to multiply by ik then transform)
751 0 : auto viewR = rho2tr_m.getView();
752 0 : const int nghostR = rho2tr_m.getNghost();
753 0 : const auto& ldomR = layoutComplex_m->getLocalNDIndex();
754 :
755 : // use temp_m as a temporary complex field
756 0 : auto view_g = temp_m.getView();
757 :
758 : // define some constants
759 0 : const scalar_type pi = Kokkos::numbers::pi_v<scalar_type>;
760 0 : const Kokkos::complex<Trhs> I = {0.0, 1.0};
761 :
762 : // define some member variables in local scope for the parallel_for
763 0 : vector_type hsize = hr_m;
764 0 : Vector<int, Dim> N = nr_m;
765 :
766 : // loop over each component (E = vector field)
767 0 : for (size_t gd = 0; gd < Dim; ++gd) {
768 : // loop over rho2tr_m to multiply by -ik (gradient in Fourier space)
769 0 : Kokkos::parallel_for(
770 0 : "Gradient - E field", rho2tr_m.getFieldRangePolicy(),
771 0 : KOKKOS_LAMBDA(const int i, const int j, const int k) {
772 : // global indices for 2N rhotr_m
773 0 : const int ig = i + ldomR[0].first() - nghostR;
774 0 : const int jg = j + ldomR[1].first() - nghostR;
775 0 : const int kg = k + ldomR[2].first() - nghostR;
776 :
777 0 : Vector<int, 3> iVec = {ig, jg, kg};
778 :
779 : scalar_type k_gd;
780 0 : const scalar_type Len = N[gd] * hsize[gd];
781 0 : const bool shift = (iVec[gd] > N[gd]);
782 0 : const bool notMid = (iVec[gd] != N[gd]);
783 :
784 0 : k_gd = notMid * (pi / Len) * (iVec[gd] - shift * 2 * N[gd]);
785 :
786 0 : view_g(i, j, k) = -(I * k_gd) * viewR(i, j, k);
787 0 : });
788 :
789 : // start a timer
790 0 : static IpplTimings::TimerRef ffte = IpplTimings::getTimer("FFT: Efield");
791 0 : IpplTimings::startTimer(ffte);
792 :
793 : // transform to get E-field
794 0 : fft_m->transform(BACKWARD, rho2_mr, temp_m);
795 :
796 0 : IpplTimings::stopTimer(ffte);
797 :
798 : // apply proper normalization
799 0 : for (unsigned int i = 0; i < Dim; ++i) {
800 0 : switch (alg) {
801 0 : case Algorithm::HOCKNEY:
802 0 : rho2_mr = rho2_mr * 2.0 * nr_m[i] * hr_m[i];
803 0 : break;
804 0 : case Algorithm::VICO:
805 : case Algorithm::BIHARMONIC:
806 0 : rho2_mr = rho2_mr * 2.0 * (1.0 / 4.0);
807 0 : break;
808 0 : case Algorithm::DCT_VICO:
809 0 : rho2_mr = rho2_mr * (1.0 / 4.0);
810 0 : break;
811 0 : default:
812 0 : throw IpplException(
813 : "FFTOpenPoissonSolver::initializeFields()",
814 : "Currently only HOCKNEY, VICO, DCT_VICO, and BIHARMONIC are "
815 : "supported for open BCs");
816 : }
817 : }
818 :
819 : // start a timer
820 0 : static IpplTimings::TimerRef edtos =
821 0 : IpplTimings::getTimer("Efield: double to phys.");
822 0 : IpplTimings::startTimer(edtos);
823 :
824 : // restrict to physical grid (N^3) and assign to LHS (E-field)
825 : // communication needed if more than one rank
826 0 : if (ranks > 1) {
827 : // COMMUNICATION
828 :
829 : // send
830 0 : const auto& lDomains1 = layout_mp->getHostLocalDomains();
831 0 : std::vector<MPI_Request> requests(0);
832 :
833 0 : for (int i = 0; i < ranks; ++i) {
834 0 : if (lDomains1[i].touches(ldom2)) {
835 0 : auto intersection = lDomains1[i].intersect(ldom2);
836 :
837 0 : solver_send(mpi::tag::OPEN_SOLVER, 0, i, intersection, ldom2, nghost2,
838 0 : view2, fd_m, requests);
839 : }
840 : }
841 :
842 : // receive
843 0 : const auto& lDomains2 = layout2_m->getHostLocalDomains();
844 :
845 0 : for (int i = 0; i < ranks; ++i) {
846 0 : if (ldom1.touches(lDomains2[i])) {
847 0 : auto intersection = ldom1.intersect(lDomains2[i]);
848 :
849 : mpi::Communicator::size_type nrecvs;
850 0 : nrecvs = intersection.size();
851 :
852 0 : buffer_type buf = Comm->getBuffer<memory_space, Trhs>(nrecvs);
853 :
854 0 : Comm->recv(i, mpi::tag::OPEN_SOLVER, fd_m, *buf, nrecvs * sizeof(Trhs),
855 : nrecvs);
856 0 : buf->resetReadPos();
857 :
858 0 : unpack(intersection, viewL, gd, fd_m, nghostL, ldom1);
859 0 : }
860 : }
861 :
862 : // wait for all messages to be received
863 0 : if (requests.size() > 0) {
864 0 : MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
865 : }
866 0 : ippl::Comm->freeAllBuffers();
867 :
868 0 : } else {
869 0 : Kokkos::parallel_for(
870 0 : "Write the E-field on physical grid", this->lhs_mp->getFieldRangePolicy(),
871 0 : KOKKOS_LAMBDA(const int i, const int j, const int k) {
872 0 : const int ig2 = i + ldom2[0].first() - nghost2;
873 0 : const int jg2 = j + ldom2[1].first() - nghost2;
874 0 : const int kg2 = k + ldom2[2].first() - nghost2;
875 :
876 0 : const int ig = i + ldom1[0].first() - nghostL;
877 0 : const int jg = j + ldom1[1].first() - nghostL;
878 0 : const int kg = k + ldom1[2].first() - nghostL;
879 :
880 : // take [0,N-1] as physical solution
881 0 : const bool isQuadrant1 = ((ig == ig2) && (jg == jg2) && (kg == kg2));
882 0 : viewL(i, j, k)[gd] = view2(i, j, k) * isQuadrant1;
883 : });
884 : }
885 0 : IpplTimings::stopTimer(edtos);
886 : }
887 0 : IpplTimings::stopTimer(efield);
888 0 : }
889 :
890 : // if user asked for Hessian, compute in Fourier domain (-k^2 multiplication)
891 0 : if (hessian) {
892 : // start a timer
893 0 : static IpplTimings::TimerRef hess = IpplTimings::getTimer("Solve: Hessian");
894 0 : IpplTimings::startTimer(hess);
895 :
896 : // get Hessian matrix view (LHS)
897 0 : auto viewH = hess_m.getView();
898 0 : const int nghostH = hess_m.getNghost();
899 :
900 : // get rho2tr_m view (as we want to multiply by -k^2 then transform)
901 0 : auto viewR = rho2tr_m.getView();
902 0 : const int nghostR = rho2tr_m.getNghost();
903 0 : const auto& ldomR = layoutComplex_m->getLocalNDIndex();
904 :
905 : // use temp_m as a temporary complex field
906 0 : auto view_g = temp_m.getView();
907 :
908 : // define some constants
909 0 : const scalar_type pi = Kokkos::numbers::pi_v<scalar_type>;
910 :
911 : // define some member variables in local scope for the parallel_for
912 0 : vector_type hsize = hr_m;
913 0 : Vector<int, Dim> N = nr_m;
914 :
915 : // loop over each component (Hessian = Matrix field)
916 0 : for (size_t row = 0; row < Dim; ++row) {
917 0 : for (size_t col = 0; col < Dim; ++col) {
918 : // loop over rho2tr_m to multiply by -k^2 (second derivative in Fourier space)
919 : // if diagonal element (row = col), do not need N/2 term = 0
920 : // else, if mixed derivative, need kVec = 0 at N/2
921 :
922 0 : Kokkos::parallel_for(
923 0 : "Hessian", rho2tr_m.getFieldRangePolicy(),
924 0 : KOKKOS_LAMBDA(const int i, const int j, const int k) {
925 : // global indices for 2N rhotr_m
926 0 : const int ig = i + ldomR[0].first() - nghostR;
927 0 : const int jg = j + ldomR[1].first() - nghostR;
928 0 : const int kg = k + ldomR[2].first() - nghostR;
929 :
930 0 : Vector<int, 3> iVec = {ig, jg, kg};
931 0 : Vector_t kVec;
932 :
933 0 : for (size_t d = 0; d < Dim; ++d) {
934 0 : const scalar_type Len = N[d] * hsize[d];
935 0 : const bool shift = (iVec[d] > N[d]);
936 0 : const bool isMid = (iVec[d] == N[d]);
937 0 : const bool notDiag = (row != col);
938 :
939 0 : kVec[d] = (1 - (notDiag * isMid)) * (pi / Len)
940 0 : * (iVec[d] - shift * 2 * N[d]);
941 : }
942 :
943 0 : view_g(i, j, k) = -(kVec[col] * kVec[row]) * viewR(i, j, k);
944 0 : });
945 :
946 : // start a timer
947 0 : static IpplTimings::TimerRef ffth = IpplTimings::getTimer("FFT: Hessian");
948 0 : IpplTimings::startTimer(ffth);
949 :
950 : // transform to get Hessian
951 0 : fft_m->transform(BACKWARD, rho2_mr, temp_m);
952 :
953 0 : IpplTimings::stopTimer(ffth);
954 :
955 : // apply proper normalization
956 0 : for (unsigned int i = 0; i < Dim; ++i) {
957 0 : switch (alg) {
958 0 : case Algorithm::HOCKNEY:
959 0 : rho2_mr = rho2_mr * 2.0 * nr_m[i] * hr_m[i];
960 0 : break;
961 0 : case Algorithm::VICO:
962 : case Algorithm::BIHARMONIC:
963 0 : rho2_mr = rho2_mr * 2.0 * (1.0 / 4.0);
964 0 : break;
965 0 : case Algorithm::DCT_VICO:
966 0 : rho2_mr = rho2_mr * (1.0 / 4.0);
967 0 : break;
968 0 : default:
969 0 : throw IpplException(
970 : "FFTOpenPoissonSolver::initializeFields()",
971 : "Currently only HOCKNEY, VICO, DCT_VICO, and BIHARMONIC are "
972 : "supported for open BCs");
973 : }
974 : }
975 :
976 : // start a timer
977 0 : static IpplTimings::TimerRef hdtos =
978 0 : IpplTimings::getTimer("Hessian: double to phys.");
979 0 : IpplTimings::startTimer(hdtos);
980 :
981 : // restrict to physical grid (N^3) and assign to Matrix field (Hessian)
982 : // communication needed if more than one rank
983 0 : if (ranks > 1) {
984 : // COMMUNICATION
985 :
986 : // send
987 0 : const auto& lDomains1 = layout_mp->getHostLocalDomains();
988 0 : std::vector<MPI_Request> requests(0);
989 :
990 0 : for (int i = 0; i < ranks; ++i) {
991 0 : if (lDomains1[i].touches(ldom2)) {
992 0 : auto intersection = lDomains1[i].intersect(ldom2);
993 :
994 0 : solver_send(mpi::tag::OPEN_SOLVER, 0, i, intersection, ldom2,
995 0 : nghost2, view2, fd_m, requests);
996 : }
997 : }
998 :
999 : // receive
1000 0 : const auto& lDomains2 = layout2_m->getHostLocalDomains();
1001 :
1002 0 : for (int i = 0; i < ranks; ++i) {
1003 0 : if (ldom1.touches(lDomains2[i])) {
1004 0 : auto intersection = ldom1.intersect(lDomains2[i]);
1005 :
1006 : mpi::Communicator::size_type nrecvs;
1007 0 : nrecvs = intersection.size();
1008 :
1009 0 : buffer_type buf = Comm->getBuffer<memory_space, Trhs>(nrecvs);
1010 :
1011 0 : Comm->recv(i, mpi::tag::OPEN_SOLVER, fd_m, *buf,
1012 : nrecvs * sizeof(Trhs), nrecvs);
1013 0 : buf->resetReadPos();
1014 :
1015 0 : unpack(intersection, viewH, fd_m, nghostH, ldom1, row, col);
1016 0 : }
1017 : }
1018 :
1019 : // wait for all messages to be received
1020 0 : if (requests.size() > 0) {
1021 0 : MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
1022 : }
1023 0 : ippl::Comm->freeAllBuffers();
1024 :
1025 0 : } else {
1026 0 : Kokkos::parallel_for(
1027 0 : "Write Hessian on physical grid", hess_m.getFieldRangePolicy(),
1028 0 : KOKKOS_LAMBDA(const int i, const int j, const int k) {
1029 0 : const int ig2 = i + ldom2[0].first() - nghost2;
1030 0 : const int jg2 = j + ldom2[1].first() - nghost2;
1031 0 : const int kg2 = k + ldom2[2].first() - nghost2;
1032 :
1033 0 : const int ig = i + ldom1[0].first() - nghostH;
1034 0 : const int jg = j + ldom1[1].first() - nghostH;
1035 0 : const int kg = k + ldom1[2].first() - nghostH;
1036 :
1037 : // take [0,N-1] as physical solution
1038 0 : const bool isQuadrant1 =
1039 0 : ((ig == ig2) && (jg == jg2) && (kg == kg2));
1040 0 : viewH(i, j, k)[row][col] = view2(i, j, k) * isQuadrant1;
1041 : });
1042 : }
1043 0 : IpplTimings::stopTimer(hdtos);
1044 : }
1045 : }
1046 0 : IpplTimings::stopTimer(hess);
1047 0 : }
1048 0 : IpplTimings::stopTimer(solve);
1049 0 : };
1050 :
1051 : ////////////////////////////////////////////////////////////////////////
1052 : // calculate FFT of the Green's function
1053 :
1054 : template <typename FieldLHS, typename FieldRHS>
1055 0 : void FFTOpenPoissonSolver<FieldLHS, FieldRHS>::greensFunction() {
1056 0 : const scalar_type pi = Kokkos::numbers::pi_v<scalar_type>;
1057 0 : grn_mr = 0.0;
1058 :
1059 0 : const int alg = this->params_m.template get<int>("algorithm");
1060 :
1061 0 : if (alg == Algorithm::VICO || alg == Algorithm::BIHARMONIC) {
1062 0 : Vector_t l(hr_m * nr_m);
1063 0 : Vector_t hs_m;
1064 0 : double L_sum(0.0);
1065 :
1066 : // compute length of the physical domain
1067 : // compute Fourier domain spacing
1068 0 : for (unsigned int i = 0; i < Dim; ++i) {
1069 0 : hs_m[i] = pi * 0.5 / l[i];
1070 0 : L_sum = L_sum + l[i] * l[i];
1071 : }
1072 :
1073 : // define the origin of the 4N grid
1074 0 : vector_type origin;
1075 :
1076 0 : for (unsigned int i = 0; i < Dim; ++i) {
1077 0 : origin[i] = -2 * nr_m[i] * pi / l[i];
1078 : }
1079 :
1080 : // set mesh for the 4N mesh
1081 0 : mesh4_m->setMeshSpacing(hs_m);
1082 :
1083 : // size of truncation window
1084 0 : L_sum = std::sqrt(L_sum);
1085 : // we choose a window 10% larger than domain (arbitrary choice)
1086 0 : L_sum = 1.1 * L_sum;
1087 :
1088 : // initialize grnL_m
1089 0 : typename CxField_gt::view_type view_g = grnL_m.getView();
1090 0 : const int nghost_g = grnL_m.getNghost();
1091 0 : const auto& ldom_g = layout4_m->getLocalNDIndex();
1092 :
1093 0 : Vector<int, Dim> size = nr_m;
1094 :
1095 : // Kokkos parallel for loop to assign analytic grnL_m
1096 0 : if (alg == Algorithm::VICO) {
1097 0 : Kokkos::parallel_for(
1098 0 : "Initialize Green's function ", grnL_m.getFieldRangePolicy(),
1099 0 : KOKKOS_LAMBDA(const int i, const int j, const int k) {
1100 : // go from local indices to global
1101 0 : const int ig = i + ldom_g[0].first() - nghost_g;
1102 0 : const int jg = j + ldom_g[1].first() - nghost_g;
1103 0 : const int kg = k + ldom_g[2].first() - nghost_g;
1104 :
1105 0 : bool isOutside = (ig > 2 * size[0] - 1);
1106 0 : const Tg t = ig * hs_m[0] + isOutside * origin[0];
1107 :
1108 0 : isOutside = (jg > 2 * size[1] - 1);
1109 0 : const Tg u = jg * hs_m[1] + isOutside * origin[1];
1110 :
1111 0 : isOutside = (kg > 2 * size[2] - 1);
1112 0 : const Tg v = kg * hs_m[2] + isOutside * origin[2];
1113 :
1114 0 : Tg s = (t * t) + (u * u) + (v * v);
1115 0 : s = Kokkos::sqrt(s);
1116 :
1117 : // assign the green's function value
1118 : // if (0,0,0), assign L^2/2 (analytical limit of sinc)
1119 :
1120 0 : const bool isOrig = ((ig == 0 && jg == 0 && kg == 0));
1121 0 : const Tg analyticLim = -L_sum * L_sum * 0.5;
1122 0 : const Tg value = -2.0 * (Kokkos::sin(0.5 * L_sum * s) / (s + isOrig * 1.0))
1123 0 : * (Kokkos::sin(0.5 * L_sum * s) / (s + isOrig * 1.0));
1124 :
1125 0 : view_g(i, j, k) = (!isOrig) * value + isOrig * analyticLim;
1126 : });
1127 0 : } else if (alg == Algorithm::BIHARMONIC) {
1128 0 : Kokkos::parallel_for(
1129 0 : "Initialize Green's function ", grnL_m.getFieldRangePolicy(),
1130 0 : KOKKOS_LAMBDA(const int i, const int j, const int k) {
1131 : // go from local indices to global
1132 0 : const int ig = i + ldom_g[0].first() - nghost_g;
1133 0 : const int jg = j + ldom_g[1].first() - nghost_g;
1134 0 : const int kg = k + ldom_g[2].first() - nghost_g;
1135 :
1136 0 : bool isOutside = (ig > 2 * size[0] - 1);
1137 0 : const Tg t = ig * hs_m[0] + isOutside * origin[0];
1138 :
1139 0 : isOutside = (jg > 2 * size[1] - 1);
1140 0 : const Tg u = jg * hs_m[1] + isOutside * origin[1];
1141 :
1142 0 : isOutside = (kg > 2 * size[2] - 1);
1143 0 : const Tg v = kg * hs_m[2] + isOutside * origin[2];
1144 :
1145 0 : Tg s = (t * t) + (u * u) + (v * v);
1146 0 : s = Kokkos::sqrt(s);
1147 :
1148 : // assign value and replace with analytic limit at origin (0,0,0)
1149 0 : const bool isOrig = ((ig == 0 && jg == 0 && kg == 0));
1150 0 : const Tg analyticLim = -L_sum * L_sum * L_sum * L_sum / 8.0;
1151 0 : const Tg value = -((2 - (L_sum * L_sum * s * s)) * Kokkos::cos(L_sum * s)
1152 0 : + 2 * L_sum * s * Kokkos::sin(L_sum * s) - 2)
1153 0 : / (2 * s * s * s * s + isOrig * 1.0);
1154 :
1155 0 : view_g(i, j, k) = (!isOrig) * value + isOrig * analyticLim;
1156 : });
1157 : }
1158 :
1159 : // start a timer
1160 0 : static IpplTimings::TimerRef fft4 = IpplTimings::getTimer("FFT: Precomputation");
1161 0 : IpplTimings::startTimer(fft4);
1162 :
1163 : // inverse Fourier transform of the green's function for precomputation
1164 0 : fft4n_m->transform(BACKWARD, grnL_m);
1165 :
1166 0 : IpplTimings::stopTimer(fft4);
1167 :
1168 : // Restrict transformed grnL_m to 2N domain after precomputation step
1169 :
1170 : // get the field data first
1171 0 : typename Field_t::view_type view = grn_mr.getView();
1172 0 : const int nghost = grn_mr.getNghost();
1173 0 : const auto& ldom = layout2_m->getLocalNDIndex();
1174 :
1175 : // start a timer
1176 0 : static IpplTimings::TimerRef ifftshift = IpplTimings::getTimer("Vico shift loop");
1177 0 : IpplTimings::startTimer(ifftshift);
1178 :
1179 : // get number of ranks to see if need communication
1180 0 : const int ranks = Comm->size();
1181 :
1182 0 : if (ranks > 1) {
1183 0 : communicateVico(size, view_g, ldom_g, nghost_g, view, ldom, nghost);
1184 : } else {
1185 : // restrict the green's function to a (2N)^3 grid from the (4N)^3 grid
1186 : using mdrange_type = Kokkos::MDRangePolicy<Kokkos::Rank<3>>;
1187 0 : Kokkos::parallel_for(
1188 : "Restrict domain of Green's function from 4N to 2N",
1189 0 : mdrange_type({nghost, nghost, nghost}, {view.extent(0) - nghost - size[0],
1190 0 : view.extent(1) - nghost - size[1],
1191 0 : view.extent(2) - nghost - size[2]}),
1192 0 : KOKKOS_LAMBDA(const int i, const int j, const int k) {
1193 : // go from local indices to global
1194 0 : const int ig = i + ldom[0].first() - nghost;
1195 0 : const int jg = j + ldom[1].first() - nghost;
1196 0 : const int kg = k + ldom[2].first() - nghost;
1197 :
1198 0 : const int ig2 = i + ldom_g[0].first() - nghost_g;
1199 0 : const int jg2 = j + ldom_g[1].first() - nghost_g;
1200 0 : const int kg2 = k + ldom_g[2].first() - nghost_g;
1201 :
1202 0 : const bool isOrig = (ig == ig2) && (jg == jg2) && (kg == kg2);
1203 0 : view(i, j, k) = isOrig * real(view_g(i, j, k));
1204 :
1205 : // Now fill the rest of the field
1206 0 : const int s = 2 * size[0] - ig - 1 - ldom_g[0].first() + nghost_g;
1207 0 : const int p = 2 * size[1] - jg - 1 - ldom_g[1].first() + nghost_g;
1208 0 : const int q = 2 * size[2] - kg - 1 - ldom_g[2].first() + nghost_g;
1209 :
1210 0 : view(s, j, k) = real(view_g(i + 1, j, k));
1211 0 : view(i, p, k) = real(view_g(i, j + 1, k));
1212 0 : view(i, j, q) = real(view_g(i, j, k + 1));
1213 0 : view(s, j, q) = real(view_g(i + 1, j, k + 1));
1214 0 : view(s, p, k) = real(view_g(i + 1, j + 1, k));
1215 0 : view(i, p, q) = real(view_g(i, j + 1, k + 1));
1216 0 : view(s, p, q) = real(view_g(i + 1, j + 1, k + 1));
1217 : });
1218 : }
1219 0 : IpplTimings::stopTimer(ifftshift);
1220 :
1221 0 : } else if (alg == Algorithm::DCT_VICO) {
1222 0 : Vector_t l(hr_m * nr_m);
1223 0 : Vector_t hs_m;
1224 0 : double L_sum(0.0);
1225 :
1226 : // compute length of the physical domain
1227 : // compute Fourier domain spacing
1228 0 : for (unsigned int i = 0; i < Dim; ++i) {
1229 0 : hs_m[i] = Kokkos::numbers::pi_v<Trhs> * 0.5 / l[i];
1230 0 : L_sum = L_sum + l[i] * l[i];
1231 : }
1232 :
1233 : // define the origin of the 2N+1 grid
1234 0 : Vector_t origin = 0.0;
1235 :
1236 : // set mesh for the 2N+1 mesh
1237 0 : mesh2n1_m->setMeshSpacing(hs_m);
1238 :
1239 : // size of truncation window
1240 0 : L_sum = Kokkos::sqrt(L_sum);
1241 0 : L_sum = 1.1 * L_sum;
1242 :
1243 : // initialize grn2n1_m
1244 0 : Vector<int, Dim> size = nr_m;
1245 :
1246 : // initialize grn2n1_m
1247 0 : typename Field_t::view_type view_g2n1 = grn2n1_m.getView();
1248 0 : const int nghost_g2n1 = grn2n1_m.getNghost();
1249 0 : const auto& ldom_g2n1 = layout2n1_m->getLocalNDIndex();
1250 :
1251 0 : Kokkos::parallel_for(
1252 0 : "Initialize 2N+1 Green's function ", grn2n1_m.getFieldRangePolicy(),
1253 0 : KOKKOS_LAMBDA(const int i, const int j, const int k) {
1254 : // go from local indices to global
1255 0 : const int ig = i + ldom_g2n1[0].first() - nghost_g2n1;
1256 0 : const int jg = j + ldom_g2n1[1].first() - nghost_g2n1;
1257 0 : const int kg = k + ldom_g2n1[2].first() - nghost_g2n1;
1258 :
1259 0 : double t = ig * hs_m[0];
1260 0 : double u = jg * hs_m[1];
1261 0 : double v = kg * hs_m[2];
1262 :
1263 0 : double s = (t * t) + (u * u) + (v * v);
1264 0 : s = Kokkos::sqrt(s);
1265 :
1266 0 : const bool isOrig = ((ig == 0 && jg == 0 && kg == 0));
1267 0 : const double val = -2.0 * (Kokkos::sin(0.5 * L_sum * s) / (s + isOrig * 1.0))
1268 0 : * (Kokkos::sin(0.5 * L_sum * s) / (s + isOrig * 1.0));
1269 0 : const double analyticLim = -L_sum * L_sum * 0.5;
1270 :
1271 0 : view_g2n1(i, j, k) = ((!isOrig) * val) + (isOrig * analyticLim);
1272 : });
1273 :
1274 : // start a timer
1275 0 : static IpplTimings::TimerRef fft4 = IpplTimings::getTimer("FFT: Precomputation");
1276 0 : IpplTimings::startTimer(fft4);
1277 :
1278 : // inverse DCT transform of 2N+1 green's function for the precomputation
1279 0 : fft2n1_m->transform(BACKWARD, grn2n1_m);
1280 :
1281 0 : IpplTimings::stopTimer(fft4);
1282 :
1283 : // Restrict transformed grn2n1_m to 2N domain after precomputation step
1284 :
1285 : // get the field data first
1286 0 : typename Field_t::view_type view = grn_mr.getView();
1287 0 : const int nghost = grn_mr.getNghost();
1288 0 : const auto& ldom = layout2_m->getLocalNDIndex();
1289 :
1290 : // start a timer
1291 0 : static IpplTimings::TimerRef ifftshift = IpplTimings::getTimer("Vico shift loop");
1292 0 : IpplTimings::startTimer(ifftshift);
1293 :
1294 : // get number of ranks to see if need communication
1295 0 : int ranks = ippl::Comm->size();
1296 :
1297 : // range type for Kokkos loop
1298 : using mdrange_type = Kokkos::MDRangePolicy<Kokkos::Rank<3>>;
1299 :
1300 0 : if (ranks > 1) {
1301 0 : communicateVico(size, view_g2n1, ldom_g2n1, nghost_g2n1, view, ldom, nghost);
1302 : } else {
1303 : // restrict the green's function to a (2N)^3 grid from the (2N+1)^3 grid
1304 0 : Kokkos::parallel_for(
1305 : "Restrict domain of Green's function from 2N+1 to 2N",
1306 0 : mdrange_type({nghost, nghost, nghost}, {view.extent(0) - nghost - size[0],
1307 0 : view.extent(1) - nghost - size[1],
1308 0 : view.extent(2) - nghost - size[2]}),
1309 0 : KOKKOS_LAMBDA(const int i, const int j, const int k) {
1310 : // go from local indices to global
1311 0 : const int ig = i + ldom[0].first() - nghost;
1312 0 : const int jg = j + ldom[1].first() - nghost;
1313 0 : const int kg = k + ldom[2].first() - nghost;
1314 :
1315 0 : const int ig2 = i + ldom_g2n1[0].first() - nghost_g2n1;
1316 0 : const int jg2 = j + ldom_g2n1[1].first() - nghost_g2n1;
1317 0 : const int kg2 = k + ldom_g2n1[2].first() - nghost_g2n1;
1318 :
1319 0 : const bool isOrig = ((ig == ig2) && (jg == jg2) && (kg == kg2));
1320 0 : view(i, j, k) = isOrig * view_g2n1(i, j, k);
1321 :
1322 : // Now fill the rest of the field
1323 0 : const int s = 2 * size[0] - ig - 1 - ldom_g2n1[0].first() + nghost_g2n1;
1324 0 : const int p = 2 * size[1] - jg - 1 - ldom_g2n1[1].first() + nghost_g2n1;
1325 0 : const int q = 2 * size[2] - kg - 1 - ldom_g2n1[2].first() + nghost_g2n1;
1326 :
1327 0 : view(s, j, k) = view_g2n1(i + 1, j, k);
1328 0 : view(i, p, k) = view_g2n1(i, j + 1, k);
1329 0 : view(i, j, q) = view_g2n1(i, j, k + 1);
1330 0 : view(s, j, q) = view_g2n1(i + 1, j, k + 1);
1331 0 : view(s, p, k) = view_g2n1(i + 1, j + 1, k);
1332 0 : view(i, p, q) = view_g2n1(i, j + 1, k + 1);
1333 0 : view(s, p, q) = view_g2n1(i + 1, j + 1, k + 1);
1334 : });
1335 : }
1336 :
1337 0 : } else {
1338 : // Hockney case
1339 :
1340 : // calculate square of the mesh spacing for each dimension
1341 0 : Vector_t hrsq(hr_m * hr_m);
1342 :
1343 : // use the grnIField_m helper field to compute Green's function
1344 0 : for (unsigned int i = 0; i < Dim; ++i) {
1345 0 : grn_mr = grn_mr + grnIField_m[i] * hrsq[i];
1346 : }
1347 :
1348 0 : grn_mr = -1.0 / (4.0 * pi * sqrt(grn_mr));
1349 :
1350 0 : typename Field_t::view_type view = grn_mr.getView();
1351 0 : const int nghost = grn_mr.getNghost();
1352 0 : const auto& ldom = layout2_m->getLocalNDIndex();
1353 :
1354 : // Kokkos parallel for loop to find (0,0,0) point and regularize
1355 0 : Kokkos::parallel_for(
1356 0 : "Regularize Green's function ", grn_mr.getFieldRangePolicy(),
1357 0 : KOKKOS_LAMBDA(const int i, const int j, const int k) {
1358 : // go from local indices to global
1359 0 : const int ig = i + ldom[0].first() - nghost;
1360 0 : const int jg = j + ldom[1].first() - nghost;
1361 0 : const int kg = k + ldom[2].first() - nghost;
1362 :
1363 : // if (0,0,0), assign to it 1/(4*pi)
1364 0 : const bool isOrig = (ig == 0 && jg == 0 && kg == 0);
1365 0 : view(i, j, k) = isOrig * (-1.0 / (4.0 * pi)) + (!isOrig) * view(i, j, k);
1366 : });
1367 0 : }
1368 :
1369 : // start a timer
1370 0 : static IpplTimings::TimerRef fftg = IpplTimings::getTimer("FFT: Green");
1371 0 : IpplTimings::startTimer(fftg);
1372 :
1373 : // perform the FFT of the Green's function for the convolution
1374 0 : fft_m->transform(FORWARD, grn_mr, grntr_m);
1375 :
1376 0 : IpplTimings::stopTimer(fftg);
1377 0 : };
1378 :
1379 : template <typename FieldLHS, typename FieldRHS>
1380 0 : void FFTOpenPoissonSolver<FieldLHS, FieldRHS>::communicateVico(
1381 : Vector<int, Dim> size, typename CxField_gt::view_type view_g,
1382 : const ippl::NDIndex<Dim> ldom_g, const int nghost_g, typename Field_t::view_type view,
1383 : const ippl::NDIndex<Dim> ldom, const int nghost) {
1384 0 : const auto& lDomains2 = layout2_m->getHostLocalDomains();
1385 0 : const auto& lDomains4 = layout4_m->getHostLocalDomains();
1386 :
1387 0 : std::vector<MPI_Request> requests(0);
1388 0 : const int ranks = Comm->size();
1389 :
1390 : // 1st step: Define 8 domains corresponding to the different quadrants
1391 0 : ippl::NDIndex<Dim> none;
1392 0 : for (unsigned i = 0; i < Dim; i++) {
1393 0 : none[i] = ippl::Index(size[i]);
1394 : }
1395 :
1396 0 : ippl::NDIndex<Dim> x;
1397 0 : x[0] = ippl::Index(size[0], 2 * size[0] - 1);
1398 0 : x[1] = ippl::Index(size[1]);
1399 0 : x[2] = ippl::Index(size[2]);
1400 :
1401 0 : ippl::NDIndex<Dim> y;
1402 0 : y[0] = ippl::Index(size[0]);
1403 0 : y[1] = ippl::Index(size[1], 2 * size[1] - 1);
1404 0 : y[2] = ippl::Index(size[2]);
1405 :
1406 0 : ippl::NDIndex<Dim> z;
1407 0 : z[0] = ippl::Index(size[0]);
1408 0 : z[1] = ippl::Index(size[1]);
1409 0 : z[2] = ippl::Index(size[2], 2 * size[2] - 1);
1410 :
1411 0 : ippl::NDIndex<Dim> xy;
1412 0 : xy[0] = ippl::Index(size[0], 2 * size[0] - 1);
1413 0 : xy[1] = ippl::Index(size[1], 2 * size[1] - 1);
1414 0 : xy[2] = ippl::Index(size[2]);
1415 :
1416 0 : ippl::NDIndex<Dim> xz;
1417 0 : xz[0] = ippl::Index(size[0], 2 * size[0] - 1);
1418 0 : xz[1] = ippl::Index(size[1]);
1419 0 : xz[2] = ippl::Index(size[2], 2 * size[2] - 1);
1420 :
1421 0 : ippl::NDIndex<Dim> yz;
1422 0 : yz[0] = ippl::Index(size[0]);
1423 0 : yz[1] = ippl::Index(size[1], 2 * size[1] - 1);
1424 0 : yz[2] = ippl::Index(size[2], 2 * size[2] - 1);
1425 :
1426 0 : ippl::NDIndex<Dim> xyz;
1427 0 : for (unsigned i = 0; i < Dim; i++) {
1428 0 : xyz[i] = ippl::Index(size[i], 2 * size[i] - 1);
1429 : }
1430 :
1431 : // 2nd step: send
1432 0 : for (int i = 0; i < ranks; ++i) {
1433 0 : auto domain2 = lDomains2[i];
1434 :
1435 0 : if (domain2.touches(none)) {
1436 0 : auto intersection = domain2.intersect(none);
1437 :
1438 0 : if (ldom_g.touches(intersection)) {
1439 0 : intersection = intersection.intersect(ldom_g);
1440 :
1441 0 : solver_send(mpi::tag::VICO_SOLVER, 0, i, intersection, ldom_g, nghost_g, view_g,
1442 0 : fd_m, requests);
1443 : }
1444 : }
1445 :
1446 0 : if (domain2.touches(x)) {
1447 0 : auto intersection = domain2.intersect(x);
1448 0 : auto xdom = ippl::Index((2 * size[0] - intersection[0].first()),
1449 0 : (2 * size[0] - intersection[0].last()), -1);
1450 :
1451 0 : ippl::NDIndex<Dim> domain4;
1452 0 : domain4[0] = xdom;
1453 0 : domain4[1] = intersection[1];
1454 0 : domain4[2] = intersection[2];
1455 :
1456 0 : if (ldom_g.touches(domain4)) {
1457 0 : intersection = ldom_g.intersect(domain4);
1458 :
1459 0 : solver_send(mpi::tag::VICO_SOLVER, 1, i, intersection, ldom_g, nghost_g, view_g,
1460 0 : fd_m, requests);
1461 : }
1462 : }
1463 :
1464 0 : if (domain2.touches(y)) {
1465 0 : auto intersection = domain2.intersect(y);
1466 0 : auto ydom = ippl::Index((2 * size[1] - intersection[1].first()),
1467 0 : (2 * size[1] - intersection[1].last()), -1);
1468 :
1469 0 : ippl::NDIndex<Dim> domain4;
1470 0 : domain4[0] = intersection[0];
1471 0 : domain4[1] = ydom;
1472 0 : domain4[2] = intersection[2];
1473 :
1474 0 : if (ldom_g.touches(domain4)) {
1475 0 : intersection = ldom_g.intersect(domain4);
1476 :
1477 0 : solver_send(mpi::tag::VICO_SOLVER, 2, i, intersection, ldom_g, nghost_g, view_g,
1478 0 : fd_m, requests);
1479 : }
1480 : }
1481 :
1482 0 : if (domain2.touches(z)) {
1483 0 : auto intersection = domain2.intersect(z);
1484 0 : auto zdom = ippl::Index((2 * size[2] - intersection[2].first()),
1485 0 : (2 * size[2] - intersection[2].last()), -1);
1486 :
1487 0 : ippl::NDIndex<Dim> domain4;
1488 0 : domain4[0] = intersection[0];
1489 0 : domain4[1] = intersection[1];
1490 0 : domain4[2] = zdom;
1491 :
1492 0 : if (ldom_g.touches(domain4)) {
1493 0 : intersection = ldom_g.intersect(domain4);
1494 :
1495 0 : solver_send(mpi::tag::VICO_SOLVER, 3, i, intersection, ldom_g, nghost_g, view_g,
1496 0 : fd_m, requests);
1497 : }
1498 : }
1499 :
1500 0 : if (domain2.touches(xy)) {
1501 0 : auto intersection = domain2.intersect(xy);
1502 0 : auto xdom = ippl::Index((2 * size[0] - intersection[0].first()),
1503 0 : (2 * size[0] - intersection[0].last()), -1);
1504 0 : auto ydom = ippl::Index((2 * size[1] - intersection[1].first()),
1505 0 : (2 * size[1] - intersection[1].last()), -1);
1506 :
1507 0 : ippl::NDIndex<Dim> domain4;
1508 0 : domain4[0] = xdom;
1509 0 : domain4[1] = ydom;
1510 0 : domain4[2] = intersection[2];
1511 :
1512 0 : if (ldom_g.touches(domain4)) {
1513 0 : intersection = ldom_g.intersect(domain4);
1514 :
1515 0 : solver_send(mpi::tag::VICO_SOLVER, 4, i, intersection, ldom_g, nghost_g, view_g,
1516 0 : fd_m, requests);
1517 : }
1518 : }
1519 :
1520 0 : if (domain2.touches(yz)) {
1521 0 : auto intersection = domain2.intersect(yz);
1522 0 : auto ydom = ippl::Index((2 * size[1] - intersection[1].first()),
1523 0 : (2 * size[1] - intersection[1].last()), -1);
1524 0 : auto zdom = ippl::Index((2 * size[2] - intersection[2].first()),
1525 0 : (2 * size[2] - intersection[2].last()), -1);
1526 :
1527 0 : ippl::NDIndex<Dim> domain4;
1528 0 : domain4[0] = intersection[0];
1529 0 : domain4[1] = ydom;
1530 0 : domain4[2] = zdom;
1531 :
1532 0 : if (ldom_g.touches(domain4)) {
1533 0 : intersection = ldom_g.intersect(domain4);
1534 0 : solver_send(mpi::tag::VICO_SOLVER, 5, i, intersection, ldom_g, nghost_g, view_g,
1535 0 : fd_m, requests);
1536 : }
1537 : }
1538 :
1539 0 : if (domain2.touches(xz)) {
1540 0 : auto intersection = domain2.intersect(xz);
1541 0 : auto xdom = ippl::Index((2 * size[0] - intersection[0].first()),
1542 0 : (2 * size[0] - intersection[0].last()), -1);
1543 0 : auto zdom = ippl::Index((2 * size[2] - intersection[2].first()),
1544 0 : (2 * size[2] - intersection[2].last()), -1);
1545 :
1546 0 : ippl::NDIndex<Dim> domain4;
1547 0 : domain4[0] = xdom;
1548 0 : domain4[1] = intersection[1];
1549 0 : domain4[2] = zdom;
1550 :
1551 0 : if (ldom_g.touches(domain4)) {
1552 0 : intersection = ldom_g.intersect(domain4);
1553 0 : solver_send(mpi::tag::VICO_SOLVER, 6, i, intersection, ldom_g, nghost_g, view_g,
1554 0 : fd_m, requests);
1555 : }
1556 : }
1557 :
1558 0 : if (domain2.touches(xyz)) {
1559 0 : auto intersection = domain2.intersect(xyz);
1560 0 : auto xdom = ippl::Index((2 * size[0] - intersection[0].first()),
1561 0 : (2 * size[0] - intersection[0].last()), -1);
1562 0 : auto ydom = ippl::Index((2 * size[1] - intersection[1].first()),
1563 0 : (2 * size[1] - intersection[1].last()), -1);
1564 0 : auto zdom = ippl::Index((2 * size[2] - intersection[2].first()),
1565 0 : (2 * size[2] - intersection[2].last()), -1);
1566 :
1567 0 : ippl::NDIndex<Dim> domain4;
1568 0 : domain4[0] = xdom;
1569 0 : domain4[1] = ydom;
1570 0 : domain4[2] = zdom;
1571 :
1572 0 : if (ldom_g.touches(domain4)) {
1573 0 : intersection = ldom_g.intersect(domain4);
1574 0 : solver_send(mpi::tag::VICO_SOLVER, 7, i, intersection, ldom_g, nghost_g, view_g,
1575 0 : fd_m, requests);
1576 : }
1577 : }
1578 : }
1579 :
1580 : // 3rd step: receive
1581 0 : for (int i = 0; i < ranks; ++i) {
1582 0 : if (ldom.touches(none)) {
1583 0 : auto intersection = ldom.intersect(none);
1584 :
1585 0 : if (lDomains4[i].touches(intersection)) {
1586 0 : intersection = intersection.intersect(lDomains4[i]);
1587 :
1588 0 : solver_recv(mpi::tag::VICO_SOLVER, 0, i, intersection, ldom, nghost, view,
1589 0 : fd_m);
1590 : }
1591 : }
1592 :
1593 0 : if (ldom.touches(x)) {
1594 0 : auto intersection = ldom.intersect(x);
1595 :
1596 0 : auto xdom = ippl::Index((2 * size[0] - intersection[0].first()),
1597 0 : (2 * size[0] - intersection[0].last()), -1);
1598 :
1599 0 : ippl::NDIndex<Dim> domain4;
1600 0 : domain4[0] = xdom;
1601 0 : domain4[1] = intersection[1];
1602 0 : domain4[2] = intersection[2];
1603 :
1604 0 : if (lDomains4[i].touches(domain4)) {
1605 0 : domain4 = lDomains4[i].intersect(domain4);
1606 0 : domain4[0] = ippl::Index(2 * size[0] - domain4[0].first(),
1607 0 : 2 * size[0] - domain4[0].last(), -1);
1608 :
1609 0 : intersection = intersection.intersect(domain4);
1610 :
1611 0 : solver_recv(mpi::tag::VICO_SOLVER, 1, i, intersection, ldom, nghost, view, fd_m,
1612 : true, false, false);
1613 : }
1614 : }
1615 :
1616 0 : if (ldom.touches(y)) {
1617 0 : auto intersection = ldom.intersect(y);
1618 :
1619 0 : auto ydom = ippl::Index((2 * size[1] - intersection[1].first()),
1620 0 : (2 * size[1] - intersection[1].last()), -1);
1621 :
1622 0 : ippl::NDIndex<Dim> domain4;
1623 0 : domain4[0] = intersection[0];
1624 0 : domain4[1] = ydom;
1625 0 : domain4[2] = intersection[2];
1626 :
1627 0 : if (lDomains4[i].touches(domain4)) {
1628 0 : domain4 = lDomains4[i].intersect(domain4);
1629 0 : domain4[1] = ippl::Index(2 * size[1] - domain4[1].first(),
1630 0 : 2 * size[1] - domain4[1].last(), -1);
1631 :
1632 0 : intersection = intersection.intersect(domain4);
1633 :
1634 0 : solver_recv(mpi::tag::VICO_SOLVER, 2, i, intersection, ldom, nghost, view, fd_m,
1635 : false, true, false);
1636 : }
1637 : }
1638 :
1639 0 : if (ldom.touches(z)) {
1640 0 : auto intersection = ldom.intersect(z);
1641 :
1642 0 : auto zdom = ippl::Index((2 * size[2] - intersection[2].first()),
1643 0 : (2 * size[2] - intersection[2].last()), -1);
1644 :
1645 0 : ippl::NDIndex<Dim> domain4;
1646 0 : domain4[0] = intersection[0];
1647 0 : domain4[1] = intersection[1];
1648 0 : domain4[2] = zdom;
1649 :
1650 0 : if (lDomains4[i].touches(domain4)) {
1651 0 : domain4 = lDomains4[i].intersect(domain4);
1652 0 : domain4[2] = ippl::Index(2 * size[2] - domain4[2].first(),
1653 0 : 2 * size[2] - domain4[2].last(), -1);
1654 :
1655 0 : intersection = intersection.intersect(domain4);
1656 :
1657 0 : solver_recv(mpi::tag::VICO_SOLVER, 3, i, intersection, ldom, nghost, view, fd_m,
1658 : false, false, true);
1659 : }
1660 : }
1661 :
1662 0 : if (ldom.touches(xy)) {
1663 0 : auto intersection = ldom.intersect(xy);
1664 :
1665 0 : auto xdom = ippl::Index((2 * size[0] - intersection[0].first()),
1666 0 : (2 * size[0] - intersection[0].last()), -1);
1667 0 : auto ydom = ippl::Index((2 * size[1] - intersection[1].first()),
1668 0 : (2 * size[1] - intersection[1].last()), -1);
1669 :
1670 0 : ippl::NDIndex<Dim> domain4;
1671 0 : domain4[0] = xdom;
1672 0 : domain4[1] = ydom;
1673 0 : domain4[2] = intersection[2];
1674 :
1675 0 : if (lDomains4[i].touches(domain4)) {
1676 0 : domain4 = lDomains4[i].intersect(domain4);
1677 0 : domain4[0] = ippl::Index(2 * size[0] - domain4[0].first(),
1678 0 : 2 * size[0] - domain4[0].last(), -1);
1679 0 : domain4[1] = ippl::Index(2 * size[1] - domain4[1].first(),
1680 0 : 2 * size[1] - domain4[1].last(), -1);
1681 :
1682 0 : intersection = intersection.intersect(domain4);
1683 :
1684 0 : solver_recv(mpi::tag::VICO_SOLVER, 4, i, intersection, ldom, nghost, view, fd_m,
1685 : true, true, false);
1686 : }
1687 : }
1688 :
1689 0 : if (ldom.touches(yz)) {
1690 0 : auto intersection = ldom.intersect(yz);
1691 :
1692 0 : auto ydom = ippl::Index((2 * size[1] - intersection[1].first()),
1693 0 : (2 * size[1] - intersection[1].last()), -1);
1694 0 : auto zdom = ippl::Index((2 * size[2] - intersection[2].first()),
1695 0 : (2 * size[2] - intersection[2].last()), -1);
1696 :
1697 0 : ippl::NDIndex<Dim> domain4;
1698 0 : domain4[0] = intersection[0];
1699 0 : domain4[1] = ydom;
1700 0 : domain4[2] = zdom;
1701 :
1702 0 : if (lDomains4[i].touches(domain4)) {
1703 0 : domain4 = lDomains4[i].intersect(domain4);
1704 0 : domain4[1] = ippl::Index(2 * size[1] - domain4[1].first(),
1705 0 : 2 * size[1] - domain4[1].last(), -1);
1706 0 : domain4[2] = ippl::Index(2 * size[2] - domain4[2].first(),
1707 0 : 2 * size[2] - domain4[2].last(), -1);
1708 :
1709 0 : intersection = intersection.intersect(domain4);
1710 :
1711 0 : solver_recv(mpi::tag::VICO_SOLVER, 5, i, intersection, ldom, nghost, view, fd_m,
1712 : false, true, true);
1713 : }
1714 : }
1715 :
1716 0 : if (ldom.touches(xz)) {
1717 0 : auto intersection = ldom.intersect(xz);
1718 :
1719 0 : auto xdom = ippl::Index((2 * size[0] - intersection[0].first()),
1720 0 : (2 * size[0] - intersection[0].last()), -1);
1721 0 : auto zdom = ippl::Index((2 * size[2] - intersection[2].first()),
1722 0 : (2 * size[2] - intersection[2].last()), -1);
1723 :
1724 0 : ippl::NDIndex<Dim> domain4;
1725 0 : domain4[0] = xdom;
1726 0 : domain4[1] = intersection[1];
1727 0 : domain4[2] = zdom;
1728 :
1729 0 : if (lDomains4[i].touches(domain4)) {
1730 0 : domain4 = lDomains4[i].intersect(domain4);
1731 0 : domain4[0] = ippl::Index(2 * size[0] - domain4[0].first(),
1732 0 : 2 * size[0] - domain4[0].last(), -1);
1733 0 : domain4[2] = ippl::Index(2 * size[2] - domain4[2].first(),
1734 0 : 2 * size[2] - domain4[2].last(), -1);
1735 :
1736 0 : intersection = intersection.intersect(domain4);
1737 :
1738 0 : solver_recv(mpi::tag::VICO_SOLVER, 6, i, intersection, ldom, nghost, view, fd_m,
1739 : true, false, true);
1740 : }
1741 : }
1742 :
1743 0 : if (ldom.touches(xyz)) {
1744 0 : auto intersection = ldom.intersect(xyz);
1745 :
1746 0 : auto xdom = ippl::Index((2 * size[0] - intersection[0].first()),
1747 0 : (2 * size[0] - intersection[0].last()), -1);
1748 0 : auto ydom = ippl::Index((2 * size[1] - intersection[1].first()),
1749 0 : (2 * size[1] - intersection[1].last()), -1);
1750 0 : auto zdom = ippl::Index((2 * size[2] - intersection[2].first()),
1751 0 : (2 * size[2] - intersection[2].last()), -1);
1752 :
1753 0 : ippl::NDIndex<Dim> domain4;
1754 0 : domain4[0] = xdom;
1755 0 : domain4[1] = ydom;
1756 0 : domain4[2] = zdom;
1757 :
1758 0 : if (lDomains4[i].touches(domain4)) {
1759 0 : domain4 = lDomains4[i].intersect(domain4);
1760 0 : domain4[0] = ippl::Index(2 * size[0] - domain4[0].first(),
1761 0 : 2 * size[0] - domain4[0].last(), -1);
1762 0 : domain4[1] = ippl::Index(2 * size[1] - domain4[1].first(),
1763 0 : 2 * size[1] - domain4[1].last(), -1);
1764 0 : domain4[2] = ippl::Index(2 * size[2] - domain4[2].first(),
1765 0 : 2 * size[2] - domain4[2].last(), -1);
1766 :
1767 0 : intersection = intersection.intersect(domain4);
1768 :
1769 0 : solver_recv(mpi::tag::VICO_SOLVER, 7, i, intersection, ldom, nghost, view, fd_m,
1770 : true, true, true);
1771 : }
1772 : }
1773 : }
1774 :
1775 0 : if (requests.size() > 0) {
1776 0 : MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
1777 : }
1778 0 : ippl::Comm->freeAllBuffers();
1779 0 : };
1780 :
1781 : // CommunicateVico for DCT_VICO (2N+1 to 2N)
1782 : template <typename FieldLHS, typename FieldRHS>
1783 0 : void FFTOpenPoissonSolver<FieldLHS, FieldRHS>::communicateVico(
1784 : Vector<int, Dim> size, typename Field_t::view_type view_g, const ippl::NDIndex<Dim> ldom_g,
1785 : const int nghost_g, typename Field_t::view_type view, const ippl::NDIndex<Dim> ldom,
1786 : const int nghost) {
1787 0 : const auto& lDomains2 = layout2_m->getHostLocalDomains();
1788 0 : const auto& lDomains2n1 = layout2n1_m->getHostLocalDomains();
1789 :
1790 0 : std::vector<MPI_Request> requests(0);
1791 0 : const int ranks = ippl::Comm->size();
1792 :
1793 : // 1st step: Define 8 domains corresponding to the different quadrants
1794 0 : ippl::NDIndex<Dim> none;
1795 0 : for (unsigned i = 0; i < Dim; i++) {
1796 0 : none[i] = ippl::Index(size[i]);
1797 : }
1798 :
1799 0 : ippl::NDIndex<Dim> x;
1800 0 : x[0] = ippl::Index(size[0], 2 * size[0] - 1);
1801 0 : x[1] = ippl::Index(size[1]);
1802 0 : x[2] = ippl::Index(size[2]);
1803 :
1804 0 : ippl::NDIndex<Dim> y;
1805 0 : y[0] = ippl::Index(size[0]);
1806 0 : y[1] = ippl::Index(size[1], 2 * size[1] - 1);
1807 0 : y[2] = ippl::Index(size[2]);
1808 :
1809 0 : ippl::NDIndex<Dim> z;
1810 0 : z[0] = ippl::Index(size[0]);
1811 0 : z[1] = ippl::Index(size[1]);
1812 0 : z[2] = ippl::Index(size[2], 2 * size[2] - 1);
1813 :
1814 0 : ippl::NDIndex<Dim> xy;
1815 0 : xy[0] = ippl::Index(size[0], 2 * size[0] - 1);
1816 0 : xy[1] = ippl::Index(size[1], 2 * size[1] - 1);
1817 0 : xy[2] = ippl::Index(size[2]);
1818 :
1819 0 : ippl::NDIndex<Dim> xz;
1820 0 : xz[0] = ippl::Index(size[0], 2 * size[0] - 1);
1821 0 : xz[1] = ippl::Index(size[1]);
1822 0 : xz[2] = ippl::Index(size[2], 2 * size[2] - 1);
1823 :
1824 0 : ippl::NDIndex<Dim> yz;
1825 0 : yz[0] = ippl::Index(size[0]);
1826 0 : yz[1] = ippl::Index(size[1], 2 * size[1] - 1);
1827 0 : yz[2] = ippl::Index(size[2], 2 * size[2] - 1);
1828 :
1829 0 : ippl::NDIndex<Dim> xyz;
1830 0 : for (unsigned i = 0; i < Dim; i++) {
1831 0 : xyz[i] = ippl::Index(size[i], 2 * size[i] - 1);
1832 : }
1833 :
1834 : // 2nd step: send
1835 0 : for (int i = 0; i < ranks; ++i) {
1836 0 : auto domain2 = lDomains2[i];
1837 :
1838 0 : if (domain2.touches(none)) {
1839 0 : auto intersection = domain2.intersect(none);
1840 :
1841 0 : if (ldom_g.touches(intersection)) {
1842 0 : intersection = intersection.intersect(ldom_g);
1843 :
1844 0 : solver_send(mpi::tag::VICO_SOLVER, 0, i, intersection, ldom_g, nghost_g, view_g,
1845 0 : fd_m, requests);
1846 : }
1847 : }
1848 :
1849 0 : if (domain2.touches(x)) {
1850 0 : auto intersection = domain2.intersect(x);
1851 0 : auto xdom = ippl::Index((2 * size[0] - intersection[0].first()),
1852 0 : (2 * size[0] - intersection[0].last()), -1);
1853 :
1854 0 : ippl::NDIndex<Dim> domain2n1;
1855 0 : domain2n1[0] = xdom;
1856 0 : domain2n1[1] = intersection[1];
1857 0 : domain2n1[2] = intersection[2];
1858 :
1859 0 : if (ldom_g.touches(domain2n1)) {
1860 0 : intersection = ldom_g.intersect(domain2n1);
1861 :
1862 0 : solver_send(mpi::tag::VICO_SOLVER, 1, i, intersection, ldom_g, nghost_g, view_g,
1863 0 : fd_m, requests);
1864 : }
1865 : }
1866 :
1867 0 : if (domain2.touches(y)) {
1868 0 : auto intersection = domain2.intersect(y);
1869 0 : auto ydom = ippl::Index((2 * size[1] - intersection[1].first()),
1870 0 : (2 * size[1] - intersection[1].last()), -1);
1871 :
1872 0 : ippl::NDIndex<Dim> domain2n1;
1873 0 : domain2n1[0] = intersection[0];
1874 0 : domain2n1[1] = ydom;
1875 0 : domain2n1[2] = intersection[2];
1876 :
1877 0 : if (ldom_g.touches(domain2n1)) {
1878 0 : intersection = ldom_g.intersect(domain2n1);
1879 :
1880 0 : solver_send(mpi::tag::VICO_SOLVER, 2, i, intersection, ldom_g, nghost_g, view_g,
1881 0 : fd_m, requests);
1882 : }
1883 : }
1884 :
1885 0 : if (domain2.touches(z)) {
1886 0 : auto intersection = domain2.intersect(z);
1887 0 : auto zdom = ippl::Index((2 * size[2] - intersection[2].first()),
1888 0 : (2 * size[2] - intersection[2].last()), -1);
1889 :
1890 0 : ippl::NDIndex<Dim> domain2n1;
1891 0 : domain2n1[0] = intersection[0];
1892 0 : domain2n1[1] = intersection[1];
1893 0 : domain2n1[2] = zdom;
1894 :
1895 0 : if (ldom_g.touches(domain2n1)) {
1896 0 : intersection = ldom_g.intersect(domain2n1);
1897 :
1898 0 : solver_send(mpi::tag::VICO_SOLVER, 3, i, intersection, ldom_g, nghost_g, view_g,
1899 0 : fd_m, requests);
1900 : }
1901 : }
1902 :
1903 0 : if (domain2.touches(xy)) {
1904 0 : auto intersection = domain2.intersect(xy);
1905 0 : auto xdom = ippl::Index((2 * size[0] - intersection[0].first()),
1906 0 : (2 * size[0] - intersection[0].last()), -1);
1907 0 : auto ydom = ippl::Index((2 * size[1] - intersection[1].first()),
1908 0 : (2 * size[1] - intersection[1].last()), -1);
1909 :
1910 0 : ippl::NDIndex<Dim> domain2n1;
1911 0 : domain2n1[0] = xdom;
1912 0 : domain2n1[1] = ydom;
1913 0 : domain2n1[2] = intersection[2];
1914 :
1915 0 : if (ldom_g.touches(domain2n1)) {
1916 0 : intersection = ldom_g.intersect(domain2n1);
1917 :
1918 0 : solver_send(mpi::tag::VICO_SOLVER, 4, i, intersection, ldom_g, nghost_g, view_g,
1919 0 : fd_m, requests);
1920 : }
1921 : }
1922 :
1923 0 : if (domain2.touches(yz)) {
1924 0 : auto intersection = domain2.intersect(yz);
1925 0 : auto ydom = ippl::Index((2 * size[1] - intersection[1].first()),
1926 0 : (2 * size[1] - intersection[1].last()), -1);
1927 0 : auto zdom = ippl::Index((2 * size[2] - intersection[2].first()),
1928 0 : (2 * size[2] - intersection[2].last()), -1);
1929 :
1930 0 : ippl::NDIndex<Dim> domain2n1;
1931 0 : domain2n1[0] = intersection[0];
1932 0 : domain2n1[1] = ydom;
1933 0 : domain2n1[2] = zdom;
1934 :
1935 0 : if (ldom_g.touches(domain2n1)) {
1936 0 : intersection = ldom_g.intersect(domain2n1);
1937 :
1938 0 : solver_send(mpi::tag::VICO_SOLVER, 5, i, intersection, ldom_g, nghost_g, view_g,
1939 0 : fd_m, requests);
1940 : }
1941 : }
1942 :
1943 0 : if (domain2.touches(xz)) {
1944 0 : auto intersection = domain2.intersect(xz);
1945 0 : auto xdom = ippl::Index((2 * size[0] - intersection[0].first()),
1946 0 : (2 * size[0] - intersection[0].last()), -1);
1947 0 : auto zdom = ippl::Index((2 * size[2] - intersection[2].first()),
1948 0 : (2 * size[2] - intersection[2].last()), -1);
1949 :
1950 0 : ippl::NDIndex<Dim> domain2n1;
1951 0 : domain2n1[0] = xdom;
1952 0 : domain2n1[1] = intersection[1];
1953 0 : domain2n1[2] = zdom;
1954 :
1955 0 : if (ldom_g.touches(domain2n1)) {
1956 0 : intersection = ldom_g.intersect(domain2n1);
1957 :
1958 0 : solver_send(mpi::tag::VICO_SOLVER, 6, i, intersection, ldom_g, nghost_g, view_g,
1959 0 : fd_m, requests);
1960 : }
1961 : }
1962 :
1963 0 : if (domain2.touches(xyz)) {
1964 0 : auto intersection = domain2.intersect(xyz);
1965 0 : auto xdom = ippl::Index((2 * size[0] - intersection[0].first()),
1966 0 : (2 * size[0] - intersection[0].last()), -1);
1967 0 : auto ydom = ippl::Index((2 * size[1] - intersection[1].first()),
1968 0 : (2 * size[1] - intersection[1].last()), -1);
1969 0 : auto zdom = ippl::Index((2 * size[2] - intersection[2].first()),
1970 0 : (2 * size[2] - intersection[2].last()), -1);
1971 :
1972 0 : ippl::NDIndex<Dim> domain2n1;
1973 0 : domain2n1[0] = xdom;
1974 0 : domain2n1[1] = ydom;
1975 0 : domain2n1[2] = zdom;
1976 :
1977 0 : if (ldom_g.touches(domain2n1)) {
1978 0 : intersection = ldom_g.intersect(domain2n1);
1979 :
1980 0 : solver_send(mpi::tag::VICO_SOLVER, 7, i, intersection, ldom_g, nghost_g, view_g,
1981 0 : fd_m, requests);
1982 : }
1983 : }
1984 : }
1985 :
1986 : // 3rd step: receive
1987 0 : for (int i = 0; i < ranks; ++i) {
1988 0 : if (ldom.touches(none)) {
1989 0 : auto intersection = ldom.intersect(none);
1990 :
1991 0 : if (lDomains2n1[i].touches(intersection)) {
1992 0 : intersection = intersection.intersect(lDomains2n1[i]);
1993 :
1994 0 : solver_recv(mpi::tag::VICO_SOLVER, 0, i, intersection, ldom, nghost, view,
1995 0 : fd_m);
1996 : }
1997 : }
1998 :
1999 0 : if (ldom.touches(x)) {
2000 0 : auto intersection = ldom.intersect(x);
2001 :
2002 0 : auto xdom = ippl::Index((2 * size[0] - intersection[0].first()),
2003 0 : (2 * size[0] - intersection[0].last()), -1);
2004 :
2005 0 : ippl::NDIndex<Dim> domain2n1;
2006 0 : domain2n1[0] = xdom;
2007 0 : domain2n1[1] = intersection[1];
2008 0 : domain2n1[2] = intersection[2];
2009 :
2010 0 : if (lDomains2n1[i].touches(domain2n1)) {
2011 0 : domain2n1 = lDomains2n1[i].intersect(domain2n1);
2012 0 : domain2n1[0] = ippl::Index(2 * size[0] - domain2n1[0].first(),
2013 0 : 2 * size[0] - domain2n1[0].last(), -1);
2014 :
2015 0 : intersection = intersection.intersect(domain2n1);
2016 :
2017 0 : solver_recv(mpi::tag::VICO_SOLVER, 1, i, intersection, ldom, nghost, view, fd_m,
2018 : true, false, false);
2019 : }
2020 : }
2021 :
2022 0 : if (ldom.touches(y)) {
2023 0 : auto intersection = ldom.intersect(y);
2024 :
2025 0 : auto ydom = ippl::Index((2 * size[1] - intersection[1].first()),
2026 0 : (2 * size[1] - intersection[1].last()), -1);
2027 :
2028 0 : ippl::NDIndex<Dim> domain2n1;
2029 0 : domain2n1[0] = intersection[0];
2030 0 : domain2n1[1] = ydom;
2031 0 : domain2n1[2] = intersection[2];
2032 :
2033 0 : if (lDomains2n1[i].touches(domain2n1)) {
2034 0 : domain2n1 = lDomains2n1[i].intersect(domain2n1);
2035 0 : domain2n1[1] = ippl::Index(2 * size[1] - domain2n1[1].first(),
2036 0 : 2 * size[1] - domain2n1[1].last(), -1);
2037 :
2038 0 : intersection = intersection.intersect(domain2n1);
2039 :
2040 0 : solver_recv(mpi::tag::VICO_SOLVER, 2, i, intersection, ldom, nghost, view, fd_m,
2041 : false, true, false);
2042 : }
2043 : }
2044 :
2045 0 : if (ldom.touches(z)) {
2046 0 : auto intersection = ldom.intersect(z);
2047 :
2048 0 : auto zdom = ippl::Index((2 * size[2] - intersection[2].first()),
2049 0 : (2 * size[2] - intersection[2].last()), -1);
2050 :
2051 0 : ippl::NDIndex<Dim> domain2n1;
2052 0 : domain2n1[0] = intersection[0];
2053 0 : domain2n1[1] = intersection[1];
2054 0 : domain2n1[2] = zdom;
2055 :
2056 0 : if (lDomains2n1[i].touches(domain2n1)) {
2057 0 : domain2n1 = lDomains2n1[i].intersect(domain2n1);
2058 0 : domain2n1[2] = ippl::Index(2 * size[2] - domain2n1[2].first(),
2059 0 : 2 * size[2] - domain2n1[2].last(), -1);
2060 :
2061 0 : intersection = intersection.intersect(domain2n1);
2062 :
2063 0 : solver_recv(mpi::tag::VICO_SOLVER, 3, i, intersection, ldom, nghost, view, fd_m,
2064 : false, false, true);
2065 : }
2066 : }
2067 :
2068 0 : if (ldom.touches(xy)) {
2069 0 : auto intersection = ldom.intersect(xy);
2070 :
2071 0 : auto xdom = ippl::Index((2 * size[0] - intersection[0].first()),
2072 0 : (2 * size[0] - intersection[0].last()), -1);
2073 0 : auto ydom = ippl::Index((2 * size[1] - intersection[1].first()),
2074 0 : (2 * size[1] - intersection[1].last()), -1);
2075 :
2076 0 : ippl::NDIndex<Dim> domain2n1;
2077 0 : domain2n1[0] = xdom;
2078 0 : domain2n1[1] = ydom;
2079 0 : domain2n1[2] = intersection[2];
2080 :
2081 0 : if (lDomains2n1[i].touches(domain2n1)) {
2082 0 : domain2n1 = lDomains2n1[i].intersect(domain2n1);
2083 0 : domain2n1[0] = ippl::Index(2 * size[0] - domain2n1[0].first(),
2084 0 : 2 * size[0] - domain2n1[0].last(), -1);
2085 0 : domain2n1[1] = ippl::Index(2 * size[1] - domain2n1[1].first(),
2086 0 : 2 * size[1] - domain2n1[1].last(), -1);
2087 :
2088 0 : intersection = intersection.intersect(domain2n1);
2089 :
2090 0 : solver_recv(mpi::tag::VICO_SOLVER, 4, i, intersection, ldom, nghost, view, fd_m,
2091 : true, true, false);
2092 : }
2093 : }
2094 :
2095 0 : if (ldom.touches(yz)) {
2096 0 : auto intersection = ldom.intersect(yz);
2097 :
2098 0 : auto ydom = ippl::Index((2 * size[1] - intersection[1].first()),
2099 0 : (2 * size[1] - intersection[1].last()), -1);
2100 0 : auto zdom = ippl::Index((2 * size[2] - intersection[2].first()),
2101 0 : (2 * size[2] - intersection[2].last()), -1);
2102 :
2103 0 : ippl::NDIndex<Dim> domain2n1;
2104 0 : domain2n1[0] = intersection[0];
2105 0 : domain2n1[1] = ydom;
2106 0 : domain2n1[2] = zdom;
2107 :
2108 0 : if (lDomains2n1[i].touches(domain2n1)) {
2109 0 : domain2n1 = lDomains2n1[i].intersect(domain2n1);
2110 0 : domain2n1[1] = ippl::Index(2 * size[1] - domain2n1[1].first(),
2111 0 : 2 * size[1] - domain2n1[1].last(), -1);
2112 0 : domain2n1[2] = ippl::Index(2 * size[2] - domain2n1[2].first(),
2113 0 : 2 * size[2] - domain2n1[2].last(), -1);
2114 :
2115 0 : intersection = intersection.intersect(domain2n1);
2116 :
2117 0 : solver_recv(mpi::tag::VICO_SOLVER, 5, i, intersection, ldom, nghost, view, fd_m,
2118 : false, true, true);
2119 : }
2120 : }
2121 :
2122 0 : if (ldom.touches(xz)) {
2123 0 : auto intersection = ldom.intersect(xz);
2124 :
2125 0 : auto xdom = ippl::Index((2 * size[0] - intersection[0].first()),
2126 0 : (2 * size[0] - intersection[0].last()), -1);
2127 0 : auto zdom = ippl::Index((2 * size[2] - intersection[2].first()),
2128 0 : (2 * size[2] - intersection[2].last()), -1);
2129 :
2130 0 : ippl::NDIndex<Dim> domain2n1;
2131 0 : domain2n1[0] = xdom;
2132 0 : domain2n1[1] = intersection[1];
2133 0 : domain2n1[2] = zdom;
2134 :
2135 0 : if (lDomains2n1[i].touches(domain2n1)) {
2136 0 : domain2n1 = lDomains2n1[i].intersect(domain2n1);
2137 0 : domain2n1[0] = ippl::Index(2 * size[0] - domain2n1[0].first(),
2138 0 : 2 * size[0] - domain2n1[0].last(), -1);
2139 0 : domain2n1[2] = ippl::Index(2 * size[2] - domain2n1[2].first(),
2140 0 : 2 * size[2] - domain2n1[2].last(), -1);
2141 :
2142 0 : intersection = intersection.intersect(domain2n1);
2143 :
2144 0 : solver_recv(mpi::tag::VICO_SOLVER, 6, i, intersection, ldom, nghost, view, fd_m,
2145 : true, false, true);
2146 : }
2147 : }
2148 :
2149 0 : if (ldom.touches(xyz)) {
2150 0 : auto intersection = ldom.intersect(xyz);
2151 :
2152 0 : auto xdom = ippl::Index((2 * size[0] - intersection[0].first()),
2153 0 : (2 * size[0] - intersection[0].last()), -1);
2154 0 : auto ydom = ippl::Index((2 * size[1] - intersection[1].first()),
2155 0 : (2 * size[1] - intersection[1].last()), -1);
2156 0 : auto zdom = ippl::Index((2 * size[2] - intersection[2].first()),
2157 0 : (2 * size[2] - intersection[2].last()), -1);
2158 :
2159 0 : ippl::NDIndex<Dim> domain2n1;
2160 0 : domain2n1[0] = xdom;
2161 0 : domain2n1[1] = ydom;
2162 0 : domain2n1[2] = zdom;
2163 :
2164 0 : if (lDomains2n1[i].touches(domain2n1)) {
2165 0 : domain2n1 = lDomains2n1[i].intersect(domain2n1);
2166 0 : domain2n1[0] = ippl::Index(2 * size[0] - domain2n1[0].first(),
2167 0 : 2 * size[0] - domain2n1[0].last(), -1);
2168 0 : domain2n1[1] = ippl::Index(2 * size[1] - domain2n1[1].first(),
2169 0 : 2 * size[1] - domain2n1[1].last(), -1);
2170 0 : domain2n1[2] = ippl::Index(2 * size[2] - domain2n1[2].first(),
2171 0 : 2 * size[2] - domain2n1[2].last(), -1);
2172 :
2173 0 : intersection = intersection.intersect(domain2n1);
2174 :
2175 0 : solver_recv(mpi::tag::VICO_SOLVER, 7, i, intersection, ldom, nghost, view, fd_m,
2176 : true, true, true);
2177 : }
2178 : }
2179 : }
2180 :
2181 0 : if (requests.size() > 0) {
2182 0 : MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
2183 : }
2184 0 : ippl::Comm->freeAllBuffers();
2185 0 : };
2186 : } // namespace ippl
|
