Block-Structured AMR Software Framework
AMReX_FFT_R2X.H
Go to the documentation of this file.
1 #ifndef AMREX_FFT_R2X_H_
2 #define AMREX_FFT_R2X_H_
3 #include <AMReX_Config.H>
4 
5 #include <AMReX_MultiFab.H>
6 #include <AMReX_FFT_Helper.H>
7 #include <algorithm>
8 #include <numeric>
9 #include <tuple>
10 
11 namespace amrex::FFT
12 {
13 
14 template <typename T> class Poisson;
15 template <typename T> class PoissonHybrid;
16 
23 template <typename T = Real>
24 class R2X
25 {
26 public:
27  using MF = std::conditional_t<std::is_same_v<T,Real>,
30 
31  template <typename U> friend class Poisson;
32  template <typename U> friend class PoissonHybrid;
33 
34  R2X (Box const& domain,
35  Array<std::pair<Boundary,Boundary>,AMREX_SPACEDIM> const& bc,
36  Info const& info = Info{});
37 
38  ~R2X ();
39 
40  R2X (R2X const&) = delete;
41  R2X (R2X &&) = delete;
42  R2X& operator= (R2X const&) = delete;
43  R2X& operator= (R2X &&) = delete;
44 
45  [[nodiscard]] T scalingFactor () const;
46 
47  template <typename F>
48  void forwardThenBackward (MF const& inmf, MF& outmf, F const& post_forward);
49 
50  // public for cuda
51  template <int dim, typename FAB, typename F>
52  void post_forward_doit (FAB* fab, F const& f);
53 
54 private:
55 
56  void forward (MF const& inmf, MF& outmf);
57  void forward (MF const& inmf, cMF& outmf);
58  void forward (MF const& inmf);
59  void backward (MF const& inmf, MF& outmf, IntVect const& ngout,
60  Periodicity const& period);
61  void backward (cMF const& inmf, MF& outmf, IntVect const& ngout,
62  Periodicity const& period);
63  void backward ();
64 
65  template <typename F>
66  void forwardThenBackward_doit (MF const& inmf, MF& outmf, F const& post_forward,
67  IntVect const& ngout = IntVect(0),
68  Periodicity const& period = Periodicity::NonPeriodic());
69 
72 
79 
80  std::unique_ptr<MultiBlockCommMetaData> m_cmd_cx2cy;
81  std::unique_ptr<MultiBlockCommMetaData> m_cmd_rx2ry;
82  std::unique_ptr<MultiBlockCommMetaData> m_cmd_cy2cz;
83  std::unique_ptr<MultiBlockCommMetaData> m_cmd_ry2rz;
84 
85  std::unique_ptr<MultiBlockCommMetaData> m_cmd_cy2cx;
86  std::unique_ptr<MultiBlockCommMetaData> m_cmd_ry2rx;
87  std::unique_ptr<MultiBlockCommMetaData> m_cmd_cz2cy;
88  std::unique_ptr<MultiBlockCommMetaData> m_cmd_rz2ry;
89 
94 
101 
102  std::unique_ptr<char,DataDeleter> m_data_1;
103  std::unique_ptr<char,DataDeleter> m_data_2;
104 
111 
113 };
114 
115 template <typename T>
116 R2X<T>::R2X (Box const& domain,
117  Array<std::pair<Boundary,Boundary>,AMREX_SPACEDIM> const& bc,
118  Info const& info)
119  : m_dom_0(domain),
120  m_bc(bc),
121  m_info(info)
122 {
123  BL_PROFILE("FFT::R2X");
124 
125  static_assert(std::is_same_v<float,T> || std::is_same_v<double,T>);
126  AMREX_ALWAYS_ASSERT(domain.smallEnd() == 0 &&
127  domain.length(0) > 1 &&
128  domain.cellCentered());
129 #if (AMREX_SPACEDIM == 3)
130  AMREX_ALWAYS_ASSERT(domain.length(1) > 1 || domain.length(2) == 1);
131 #else
133 #endif
134  for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
135  AMREX_ALWAYS_ASSERT(domain.length(idim) > 1);
136  if (bc[idim].first == Boundary::periodic ||
137  bc[idim].second == Boundary::periodic) {
138  AMREX_ALWAYS_ASSERT(bc[idim].first == bc[idim].second);
139  }
140  }
141 
142  int myproc = ParallelContext::MyProcSub();
144 
145  //
146  // make data containers
147  //
148 
149  m_dom_rx = m_dom_0;
150  auto bax = amrex::decompose(m_dom_rx, nprocs, {AMREX_D_DECL(false,true,true)});
152  m_rx.define(bax, dmx, 1, 0, MFInfo().SetAlloc(false));
153 
154  // x-direction
155  if (bc[0].first == Boundary::periodic) {
156  // x-fft: r2c(m_rx->m_cx)
157  m_dom_cx = Box(IntVect(0), IntVect(AMREX_D_DECL(domain.length(0)/2,
158  domain.bigEnd(1),
159  domain.bigEnd(2))));
160  BoxList bl = bax.boxList();
161  for (auto & b : bl) {
162  b.setBig(0, m_dom_cx.bigEnd(0));
163  }
164  BoxArray cbax(std::move(bl));
165  m_cx.define(cbax, dmx, 1, 0, MFInfo().SetAlloc(false));
166  } // else: x-fft: r2r(m_rx)
167 
168 #if (AMREX_SPACEDIM >= 2)
169  if (domain.length(1) > 1) {
170  if (! m_cx.empty()) {
171  // copy(m_cx->m_cy)
173  m_dom_cx.bigEnd(0),
174  m_dom_cx.bigEnd(2))));
175  auto ba = amrex::decompose(m_dom_cy, nprocs, {AMREX_D_DECL(false,true,true)});
177  if (ba.size() == m_cx.size()) {
178  dm = m_cx.DistributionMap();
179  } else {
180  dm = detail::make_iota_distromap(ba.size());
181  }
182  m_cy.define(ba, dm, 1, 0, MFInfo().SetAlloc(false));
183  // if bc[1] is periodic:
184  // c2c(m_cy->m_cy)
185  // else:
186  // r2r(m_cy.re) & r2r(m_cy.im)
187  } else {
188  // copy(m_rx->m_ry)
190  m_dom_rx.bigEnd(0),
191  m_dom_rx.bigEnd(2))));
192  auto ba = amrex::decompose(m_dom_ry, nprocs, {AMREX_D_DECL(false,true,true)});
194  if (ba.size() == m_rx.size()) {
195  dm = m_rx.DistributionMap();
196  } else {
197  dm = detail::make_iota_distromap(ba.size());
198  }
199  m_ry.define(ba, dm, 1, 0, MFInfo().SetAlloc(false));
200  // if bc[1] is periodic:
201  // r2c(m_ry->m_cy)
202  // else:
203  // r2r(m_ry)
204  if (bc[1].first == Boundary::periodic) {
206  m_dom_ry.bigEnd(1),
207  m_dom_ry.bigEnd(2))));
208  BoxList bl = ba.boxList();
209  for (auto & b : bl) {
210  b.setBig(0, m_dom_cy.bigEnd(0));
211  }
212  BoxArray cba(std::move(bl));
213  m_cy.define(cba, dm, 1, 0, MFInfo().SetAlloc(false));
214  }
215  }
216  }
217 #endif
218 
219 #if (AMREX_SPACEDIM == 3)
220  if (domain.length(2) > 1 && !m_info.batch_mode) {
221  if (! m_cy.empty()) {
222  // copy(m_cy, m_cz)
224  m_dom_cy.bigEnd(1),
225  m_dom_cy.bigEnd(0))));
226  auto ba = amrex::decompose(m_dom_cz, nprocs, {AMREX_D_DECL(false,true,true)});
228  if (ba.size() == m_cy.size()) {
229  dm = m_cy.DistributionMap();
230  } else {
231  dm = detail::make_iota_distromap(ba.size());
232  }
233  m_cz.define(ba, dm, 1, 0, MFInfo().SetAlloc(false));
234  // if bc[2] is periodic:
235  // c2c(m_cz->m_cz)
236  // else:
237  // r2r(m_cz.re) & r2r(m_cz.im)
238  } else {
239  // copy(m_ry, m_rz)
241  m_dom_ry.bigEnd(1),
242  m_dom_ry.bigEnd(0))));
243  auto ba = amrex::decompose(m_dom_rz, nprocs, {AMREX_D_DECL(false,true,true)});
245  if (ba.size() == m_ry.size()) {
246  dm = m_ry.DistributionMap();
247  } else {
248  dm = detail::make_iota_distromap(ba.size());
249  }
250  m_rz.define(ba, dm, 1, 0, MFInfo().SetAlloc(false));
251  // if bc[2] is periodic:
252  // r2c(m_rz->m_cz)
253  // else:
254  // r2r(m_rz)
255  if (bc[2].first == Boundary::periodic) {
257  m_dom_rz.bigEnd(1),
258  m_dom_rz.bigEnd(2))));
259  BoxList bl = ba.boxList();
260  for (auto & b : bl) {
261  b.setBig(0, m_dom_cz.bigEnd(0));
262  }
263  BoxArray cba(std::move(bl));
264  m_cz.define(cba, dm, 1, 0, MFInfo().SetAlloc(false));
265  }
266  }
267  }
268 #endif
269 
270  // There are several different execution paths.
271  //
272  // (1) x-r2c(m_rx->m_cx), copy(m_cx->m_cy), y-fft(m_cy),
273  // copy(m_cy->m_cz), z-fft(m_cz)
274  // In this case, we have m_rx, m_cx, m_cy, & m_cz.
275  // we can alias(m_rx,m_cy) and alias(m_cx,m_cz).
276  //
277  // (2) x_r2r(m_rx), copy(m_rx->m_ry), y-r2c(m_ry->m_cy),
278  // copy(m_cy->m_cz), z-fft(m_cz)
279  // In this case, we have m_rx, m_ry, m_cy, & m_cz.
280  // We can alias(m_rx,m_cy) and alias(m_ry,m_cz).
281  //
282  // (3) x_r2r(m_rx), copy(m_rx->m_ry), y-r2r(m_ry),
283  // copy(m_ry->m_rz), z-r2c(m_rz->m_rz)
284  // In this case, we have m_rx, m_ry, m_rz, & m_cz
285  // We can alias(m_rx,m_rz) and alias(m_ry,m_cz)
286  //
287  // (4) x_r2r(m_rx), copy(m_rx->m_ry), y-r2r(m_ry),
288  // copy(m_ry->m_rz), z-r2r(m_rz)
289  // In this case, we have m_rx, m_ry, & m_rz.
290  // We can alias(m_rx,m_rz).
291 
292  if (! m_cx.empty()) {
295  } else if (! m_cy.empty()) {
298  } else if (! m_cz.empty()) {
301  } else {
303  m_data_2 = detail::make_mfs_share(m_ry, m_cz); // It's okay m_cz is empty.
304  }
305 
306  //
307  // make copiers
308  //
309 
310 #if (AMREX_SPACEDIM >= 2)
311  if (domain.length(1) > 1) {
312  if (! m_cx.empty()) {
313  // copy(m_cx->m_cy)
314  m_cmd_cx2cy = std::make_unique<MultiBlockCommMetaData>
316  m_cmd_cy2cx = std::make_unique<MultiBlockCommMetaData>
318  } else {
319  // copy(m_rx->m_ry)
320  m_cmd_rx2ry = std::make_unique<MultiBlockCommMetaData>
322  m_cmd_ry2rx = std::make_unique<MultiBlockCommMetaData>
324  }
325  }
326 #endif
327 
328 #if (AMREX_SPACEDIM == 3)
329  if (domain.length(2) > 1 && !m_info.batch_mode) {
330  if (! m_cy.empty()) {
331  // copy(m_cy, m_cz)
332  m_cmd_cy2cz = std::make_unique<MultiBlockCommMetaData>
334  m_cmd_cz2cy = std::make_unique<MultiBlockCommMetaData>
336  } else {
337  // copy(m_ry, m_rz)
338  m_cmd_ry2rz = std::make_unique<MultiBlockCommMetaData>
340  m_cmd_rz2ry = std::make_unique<MultiBlockCommMetaData>
342  }
343  }
344 #endif
345 
346  //
347  // make plans
348  //
349 
350  using VendorComplex = typename Plan<T>::VendorComplex;
351 
352  if (myproc < m_rx.size())
353  {
354  Box const& box = m_rx.box(myproc);
355  auto* pf = m_rx[myproc].dataPtr();
356  if (bc[0].first == Boundary::periodic) {
357  auto* pb = (VendorComplex*) m_cx[myproc].dataPtr();
358  m_fft_fwd_x.template init_r2c<Direction::forward>(box, pf, pb);
359 #if defined(AMREX_USE_SYCL)
361 #else
362  m_fft_bwd_x.template init_r2c<Direction::backward>(box, pf, pb);
363 #endif
364  } else {
365  m_fft_fwd_x.template init_r2r<Direction::forward>(box, pf, bc[0]);
366 #if defined(AMREX_USE_GPU)
367  if ((bc[0].first == Boundary::even && bc[0].second == Boundary::odd) ||
368  (bc[0].first == Boundary::odd && bc[0].second == Boundary::even)) {
370  } else
371 #endif
372  {
373  m_fft_bwd_x.template init_r2r<Direction::backward>(box, pf, bc[0]);
374  }
375  }
376  }
377 
378 #if (AMREX_SPACEDIM >= 2)
379  if (m_ry.empty() && m_bc[1].first == Boundary::periodic) {
380  if (myproc < m_cy.size()) {
381  Box const& box = m_cy.box(myproc);
382  auto* p = (VendorComplex *)m_cy[myproc].dataPtr();
383  m_fft_fwd_y.template init_c2c<Direction::forward>(box, p);
384 #if defined(AMREX_USE_SYCL)
386 #else
387  m_fft_bwd_y.template init_c2c<Direction::backward>(box, p);
388 #endif
389  }
390  } else if (!m_ry.empty() && m_bc[1].first == Boundary::periodic) {
391  if (myproc < m_ry.size()) {
392  Box const& box = m_ry.box(myproc);
393  auto* pr = m_ry[myproc].dataPtr();
394  auto* pc = (VendorComplex*)m_cy[myproc].dataPtr();
395  m_fft_fwd_y.template init_r2c<Direction::forward>(box, pr, pc);
396 #if defined(AMREX_USE_SYCL)
398 #else
399  m_fft_bwd_y.template init_r2c<Direction::backward>(box, pr, pc);
400 #endif
401  }
402  } else if (!m_cy.empty()) {
403  if (myproc < m_cy.size()) {
404  Box const& box = m_cy.box(myproc);
405  auto* p = (VendorComplex*) m_cy[myproc].dataPtr();
406  m_fft_fwd_y.template init_r2r<Direction::forward>(box, p, bc[1]);
407 #if defined(AMREX_USE_GPU)
408  if ((bc[1].first == Boundary::even && bc[1].second == Boundary::odd) ||
409  (bc[1].first == Boundary::odd && bc[1].second == Boundary::even)) {
411  } else
412 #endif
413  {
414  m_fft_bwd_y.template init_r2r<Direction::backward>(box, p, bc[1]);
415  }
416  }
417  } else {
418  if (myproc < m_ry.size()) {
419  Box const& box = m_ry.box(myproc);
420  auto* p = m_ry[myproc].dataPtr();
421  m_fft_fwd_y.template init_r2r<Direction::forward>(box, p, bc[1]);
422 #if defined(AMREX_USE_GPU)
423  if ((bc[1].first == Boundary::even && bc[1].second == Boundary::odd) ||
424  (bc[1].first == Boundary::odd && bc[1].second == Boundary::even)) {
426  } else
427 #endif
428  {
429  m_fft_bwd_y.template init_r2r<Direction::backward>(box, p, bc[1]);
430  }
431  }
432  }
433 #endif
434 
435 #if (AMREX_SPACEDIM == 3)
436  if (m_rz.empty() && m_bc[2].first == Boundary::periodic) {
437  if (myproc < m_cz.size()) {
438  Box const& box = m_cz.box(myproc);
439  auto* p = (VendorComplex*)m_cz[myproc].dataPtr();
440  m_fft_fwd_z.template init_c2c<Direction::forward>(box, p);
441 #if defined(AMREX_USE_SYCL)
443 #else
444  m_fft_bwd_z.template init_c2c<Direction::backward>(box, p);
445 #endif
446  }
447  } else if (!m_rz.empty() && m_bc[2].first == Boundary::periodic) {
448  if (myproc < m_rz.size()) {
449  Box const& box = m_rz.box(myproc);
450  auto* pr = m_rz[myproc].dataPtr();
451  auto* pc = (VendorComplex*)m_cz[myproc].dataPtr();
452  m_fft_fwd_z.template init_r2c<Direction::forward>(box, pr, pc);
453 #if defined(AMREX_USE_SYCL)
455 #else
456  m_fft_bwd_z.template init_r2c<Direction::backward>(box, pr, pc);
457 #endif
458  }
459  } else if (!m_cz.empty()) {
460  if (myproc < m_cz.size()) {
461  Box const& box = m_cz.box(myproc);
462  auto* p = (VendorComplex*) m_cz[myproc].dataPtr();
463  m_fft_fwd_z.template init_r2r<Direction::forward>(box, p, bc[2]);
464 #if defined(AMREX_USE_GPU)
465  if ((bc[2].first == Boundary::even && bc[2].second == Boundary::odd) ||
466  (bc[2].first == Boundary::odd && bc[2].second == Boundary::even)) {
468  } else
469 #endif
470  {
471  m_fft_bwd_z.template init_r2r<Direction::backward>(box, p, bc[2]);
472  }
473  }
474  } else {
475  if (myproc < m_rz.size()) {
476  Box const& box = m_rz.box(myproc);
477  auto* p = m_rz[myproc].dataPtr();
478  m_fft_fwd_z.template init_r2r<Direction::forward>(box, p, bc[2]);
479 #if defined(AMREX_USE_GPU)
480  if ((bc[2].first == Boundary::even && bc[2].second == Boundary::odd) ||
481  (bc[2].first == Boundary::odd && bc[2].second == Boundary::even)) {
483  } else
484 #endif
485  {
486  m_fft_bwd_z.template init_r2r<Direction::backward>(box, p, bc[2]);
487  }
488  }
489  }
490 #endif
491 }
492 
493 template <typename T>
495 {
496  if (m_fft_bwd_x.plan != m_fft_fwd_x.plan) {
497  m_fft_bwd_x.destroy();
498  }
499  if (m_fft_bwd_y.plan != m_fft_fwd_y.plan) {
500  m_fft_bwd_y.destroy();
501  }
502  if (m_fft_bwd_z.plan != m_fft_fwd_z.plan) {
503  m_fft_bwd_z.destroy();
504  }
505  m_fft_fwd_x.destroy();
506  m_fft_fwd_y.destroy();
507  m_fft_fwd_z.destroy();
508 }
509 
510 template <typename T>
512 {
513  Long r = 1;
514  int ndims = m_info.batch_mode ? AMREX_SPACEDIM-1 : AMREX_SPACEDIM;
515  for (int idim = 0; idim < ndims; ++idim) {
516  r *= m_dom_0.length(idim);
517  if (m_bc[idim].first != Boundary::periodic && (m_dom_0.length(idim) > 1)) {
518  r *= 2;
519  }
520  }
521  return T(1)/T(r);
522 }
523 
524 template <typename T>
525 template <typename F>
526 void R2X<T>::forwardThenBackward (MF const& inmf, MF& outmf, F const& post_forward)
527 {
528  forwardThenBackward_doit(inmf, outmf, post_forward);
529 }
530 
531 template <typename T>
532 template <typename F>
533 void R2X<T>::forwardThenBackward_doit (MF const& inmf, MF& outmf,
534  F const& post_forward,
535  IntVect const& ngout,
536  Periodicity const& period)
537 {
538  BL_PROFILE("FFT::R2X::forwardbackward");
539 
540  this->forward(inmf);
541 
542  // post-forward
543 
544  int actual_dim = AMREX_SPACEDIM;
545 #if (AMREX_SPACEDIM >= 2)
546  if (m_dom_0.length(1) == 1) { actual_dim = 1; }
547 #endif
548 #if (AMREX_SPACEDIM == 3)
549  if ((m_dom_0.length(2) == 1) && (m_dom_0.length(1) > 1)) { actual_dim = 2; }
550 #endif
551 
552  if (actual_dim == 1) {
553  if (m_cx.empty()) {
554  post_forward_doit<0>(detail::get_fab(m_rx), post_forward);
555  } else {
556  post_forward_doit<0>(detail::get_fab(m_cx), post_forward);
557  }
558  }
559 #if (AMREX_SPACEDIM >= 2)
560  else if (actual_dim == 2) {
561  if (m_cy.empty()) {
562  post_forward_doit<1>(detail::get_fab(m_ry), post_forward);
563  } else {
564  post_forward_doit<1>(detail::get_fab(m_cy), post_forward);
565  }
566  }
567 #endif
568 #if (AMREX_SPACEDIM == 3)
569  else if (actual_dim == 3) {
570  if (m_cz.empty()) {
571  post_forward_doit<2>(detail::get_fab(m_rz), post_forward);
572  } else {
573  post_forward_doit<2>(detail::get_fab(m_cz), post_forward);
574  }
575  }
576 #endif
577 
578  this->backward();
579 
580  outmf.ParallelCopy(m_rx, 0, 0, 1, IntVect(0), ngout, period);
581 }
582 
583 template <typename T>
584 void R2X<T>::forward (MF const& inmf)
585 {
586  BL_PROFILE("FFT::R2X::forward");
587 
588  m_rx.ParallelCopy(inmf, 0, 0, 1);
589  if (m_bc[0].first == Boundary::periodic) {
590  m_fft_fwd_x.template compute_r2c<Direction::forward>();
591  } else {
592  m_fft_fwd_x.template compute_r2r<Direction::forward>();
593  }
594 
595 #if (AMREX_SPACEDIM >= 2)
596  if ( m_cmd_cx2cy) {
597  ParallelCopy(m_cy, m_cx, *m_cmd_cx2cy, 0, 0, 1, m_dtos_x2y);
598  } else if ( m_cmd_rx2ry) {
599  ParallelCopy(m_ry, m_rx, *m_cmd_rx2ry, 0, 0, 1, m_dtos_x2y);
600  }
601  if (m_bc[1].first != Boundary::periodic)
602  {
603  m_fft_fwd_y.template compute_r2r<Direction::forward>();
604  }
605  else if (m_bc[0].first == Boundary::periodic)
606  {
607  m_fft_fwd_y.template compute_c2c<Direction::forward>();
608  }
609  else
610  {
611  m_fft_fwd_y.template compute_r2c<Direction::forward>();
612  }
613 #endif
614 
615 #if (AMREX_SPACEDIM == 3)
616  if ( m_cmd_cy2cz) {
617  ParallelCopy(m_cz, m_cy, *m_cmd_cy2cz, 0, 0, 1, m_dtos_y2z);
618  } else if ( m_cmd_ry2rz) {
619  ParallelCopy(m_rz, m_ry, *m_cmd_ry2rz, 0, 0, 1, m_dtos_y2z);
620  }
621  if (m_bc[2].first != Boundary::periodic)
622  {
623  m_fft_fwd_z.template compute_r2r<Direction::forward>();
624  }
625  else if (m_bc[0].first == Boundary::periodic ||
626  m_bc[1].first == Boundary::periodic)
627  {
628  m_fft_fwd_z.template compute_c2c<Direction::forward>();
629  }
630  else
631  {
632  m_fft_fwd_z.template compute_r2c<Direction::forward>();
633  }
634 #endif
635 }
636 
637 template <typename T>
638 void R2X<T>::forward (MF const& inmf, MF& outmf)
639 {
640  this->forward(inmf);
641 
642 #if (AMREX_SPACEDIM == 3)
643  if (m_info.batch_mode) {
644  if (m_cy.empty()) {
645  ParallelCopy(outmf, m_dom_rx, m_ry, 0, 0, 1, IntVect(0), Swap01{});
646  } else {
647  amrex::Abort("R2X::forward(MF,MF): How did this happen?");
648  }
649  } else
650 #endif
651  {
652  amrex::ignore_unused(outmf);
653  amrex::Abort("R2X::forward(MF,MF): TODO");
654  }
655 }
656 
657 template <typename T>
658 void R2X<T>::forward (MF const& inmf, cMF& outmf)
659 {
660  this->forward(inmf);
661 
662 #if (AMREX_SPACEDIM == 3)
663  if (m_info.batch_mode) {
664  if (!m_cy.empty()) {
665  auto lo = m_dom_cy.smallEnd();
666  auto hi = m_dom_cy.bigEnd();
667  std::swap(lo[0],lo[1]);
668  std::swap(hi[0],hi[1]);
669  Box dom(lo,hi);
670  ParallelCopy(outmf, dom, m_cy, 0, 0, 1, IntVect(0), Swap01{});
671  } else {
672  amrex::Abort("R2X::forward(MF,cMF): How did this happen?");
673  }
674  } else
675 #endif
676  {
677  amrex::ignore_unused(outmf);
678  amrex::Abort("R2X::forward(MF,cMF): TODO");
679  }
680 }
681 
682 template <typename T>
684 {
685  BL_PROFILE("FFT::R2X::backward");
686 
687 #if (AMREX_SPACEDIM == 3)
688  if (m_bc[2].first != Boundary::periodic)
689  {
690  m_fft_bwd_z.template compute_r2r<Direction::backward>();
691  }
692  else if (m_bc[0].first == Boundary::periodic ||
693  m_bc[1].first == Boundary::periodic)
694  {
695  m_fft_bwd_z.template compute_c2c<Direction::backward>();
696  }
697  else
698  {
699  m_fft_bwd_z.template compute_r2c<Direction::backward>();
700  }
701  if ( m_cmd_cz2cy) {
702  ParallelCopy(m_cy, m_cz, *m_cmd_cz2cy, 0, 0, 1, m_dtos_z2y);
703  } else if ( m_cmd_rz2ry) {
704  ParallelCopy(m_ry, m_rz, *m_cmd_rz2ry, 0, 0, 1, m_dtos_z2y);
705  }
706 #endif
707 
708 #if (AMREX_SPACEDIM >= 2)
709  if (m_bc[1].first != Boundary::periodic)
710  {
711  m_fft_bwd_y.template compute_r2r<Direction::backward>();
712  }
713  else if (m_bc[0].first == Boundary::periodic)
714  {
715  m_fft_bwd_y.template compute_c2c<Direction::backward>();
716  }
717  else
718  {
719  m_fft_bwd_y.template compute_r2c<Direction::backward>();
720  }
721  if ( m_cmd_cy2cx) {
722  ParallelCopy(m_cx, m_cy, *m_cmd_cy2cx, 0, 0, 1, m_dtos_y2x);
723  } else if ( m_cmd_ry2rx) {
724  ParallelCopy(m_rx, m_ry, *m_cmd_ry2rx, 0, 0, 1, m_dtos_y2x);
725  }
726 #endif
727 
728  if (m_bc[0].first == Boundary::periodic) {
729  m_fft_bwd_x.template compute_r2c<Direction::backward>();
730  } else {
731  m_fft_bwd_x.template compute_r2r<Direction::backward>();
732  }
733 }
734 
735 template <typename T>
736 void R2X<T>::backward (MF const& inmf, MF& outmf, IntVect const& ngout,
737  Periodicity const& period)
738 {
739 #if (AMREX_SPACEDIM == 3)
740  if (m_info.batch_mode) {
741  if (m_cy.empty()) {
742  ParallelCopy(m_ry, m_dom_ry, inmf, 0, 0, 1, IntVect(0), Swap01{});
743  } else {
744  amrex::Abort("R2X::backward(MF,MF): How did this happen?");
745  }
746  } else
747 #endif
748  {
749  amrex::ignore_unused(inmf,outmf,ngout,period);
750  amrex::Abort("R2X::backward(MF,MF): TODO");
751  }
752 
753  this->backward();
754 
755  outmf.ParallelCopy(m_rx, 0, 0, 1, IntVect(0), ngout, period);
756 }
757 
758 template <typename T>
759 void R2X<T>::backward (cMF const& inmf, MF& outmf, IntVect const& ngout,
760  Periodicity const& period)
761 {
762 #if (AMREX_SPACEDIM == 3)
763  if (m_info.batch_mode) {
764  if (!m_cy.empty()) {
765  ParallelCopy(m_cy, m_dom_cy, inmf, 0, 0, 1, IntVect(0), Swap01{});
766  } else {
767  amrex::Abort("R2X::backward(cMF,MF): How did this happen?");
768  }
769  } else
770 #endif
771  {
772  amrex::ignore_unused(inmf,outmf,ngout,period);
773  amrex::Abort("R2X::backward(cMF,MF): TODO");
774  }
775 
776  this->backward();
777 
778  outmf.ParallelCopy(m_rx, 0, 0, 1, IntVect(0), ngout, period);
779 }
780 
781 template <typename T>
782 template <int dim, typename FAB, typename F>
783 void R2X<T>::post_forward_doit (FAB* fab, F const& f)
784 {
785  if (fab) {
786  auto const& a = fab->array();
787  ParallelFor(fab->box(),
788  [f=f,a=a] AMREX_GPU_DEVICE (int i, int j, int k)
789  {
790  if constexpr (dim == 0) {
791  f(i,j,k,a(i,j,k));
792  } else if constexpr (dim == 1) {
793  f(j,i,k,a(i,j,k));
794  } else {
795  f(j,k,i,a(i,j,k));
796  }
797  });
798  }
799 }
800 
801 }
802 
803 #endif
#define BL_PROFILE(a)
Definition: AMReX_BLProfiler.H:551
#define AMREX_ALWAYS_ASSERT(EX)
Definition: AMReX_BLassert.H:50
#define AMREX_GPU_DEVICE
Definition: AMReX_GpuQualifiers.H:18
#define AMREX_D_DECL(a, b, c)
Definition: AMReX_SPACE.H:104
A collection of Boxes stored in an Array.
Definition: AMReX_BoxArray.H:550
A class for managing a List of Boxes that share a common IndexType. This class implements operations ...
Definition: AMReX_BoxList.H:52
AMREX_GPU_HOST_DEVICE const IntVectND< dim > & smallEnd() const &noexcept
Get the smallend of the BoxND.
Definition: AMReX_Box.H:105
AMREX_GPU_HOST_DEVICE const IntVectND< dim > & bigEnd() const &noexcept
Get the bigend.
Definition: AMReX_Box.H:116
AMREX_GPU_HOST_DEVICE IntVectND< dim > length() const noexcept
Return the length of the BoxND.
Definition: AMReX_Box.H:146
AMREX_GPU_HOST_DEVICE bool cellCentered() const noexcept
Returns true if BoxND is cell-centered in all indexing directions.
Definition: AMReX_Box.H:319
Calculates the distribution of FABs to MPI processes.
Definition: AMReX_DistributionMapping.H:41
3D Poisson solver for periodic, Dirichlet & Neumann boundaries in the first two dimensions,...
Definition: AMReX_FFT_Poisson.H:106
Poisson solver for periodic, Dirichlet & Neumann boundaries using FFT.
Definition: AMReX_FFT_Poisson.H:22
Discrete Fourier Transform.
Definition: AMReX_FFT_R2X.H:25
MF m_rz
Definition: AMReX_FFT_R2X.H:97
Array< std::pair< Boundary, Boundary >, AMREX_SPACEDIM > m_bc
Definition: AMReX_FFT_R2X.H:71
void post_forward_doit(FAB *fab, F const &f)
Definition: AMReX_FFT_R2X.H:783
Swap01 m_dtos_x2y
Definition: AMReX_FFT_R2X.H:90
std::conditional_t< std::is_same_v< T, Real >, MultiFab, FabArray< BaseFab< T > > > MF
Definition: AMReX_FFT_R2X.H:28
std::unique_ptr< MultiBlockCommMetaData > m_cmd_rz2ry
Definition: AMReX_FFT_R2X.H:88
MF m_ry
Definition: AMReX_FFT_R2X.H:96
Box m_dom_cz
Definition: AMReX_FFT_R2X.H:110
Swap02 m_dtos_y2z
Definition: AMReX_FFT_R2X.H:92
~R2X()
Definition: AMReX_FFT_R2X.H:494
void forward(MF const &inmf, MF &outmf)
Definition: AMReX_FFT_R2X.H:638
void forwardThenBackward(MF const &inmf, MF &outmf, F const &post_forward)
Definition: AMReX_FFT_R2X.H:526
std::unique_ptr< MultiBlockCommMetaData > m_cmd_ry2rx
Definition: AMReX_FFT_R2X.H:86
Info m_info
Definition: AMReX_FFT_R2X.H:112
std::unique_ptr< MultiBlockCommMetaData > m_cmd_cz2cy
Definition: AMReX_FFT_R2X.H:87
std::unique_ptr< MultiBlockCommMetaData > m_cmd_cy2cz
Definition: AMReX_FFT_R2X.H:82
Plan< T > m_fft_fwd_y
Definition: AMReX_FFT_R2X.H:75
MF m_rx
Definition: AMReX_FFT_R2X.H:95
std::unique_ptr< MultiBlockCommMetaData > m_cmd_ry2rz
Definition: AMReX_FFT_R2X.H:83
Box m_dom_rx
Definition: AMReX_FFT_R2X.H:105
std::unique_ptr< MultiBlockCommMetaData > m_cmd_rx2ry
Definition: AMReX_FFT_R2X.H:81
Swap01 m_dtos_y2x
Definition: AMReX_FFT_R2X.H:91
Box m_dom_cy
Definition: AMReX_FFT_R2X.H:109
R2X(R2X &&)=delete
T scalingFactor() const
Definition: AMReX_FFT_R2X.H:511
Box m_dom_ry
Definition: AMReX_FFT_R2X.H:106
void backward()
Definition: AMReX_FFT_R2X.H:683
Plan< T > m_fft_bwd_z
Definition: AMReX_FFT_R2X.H:78
void forwardThenBackward_doit(MF const &inmf, MF &outmf, F const &post_forward, IntVect const &ngout=IntVect(0), Periodicity const &period=Periodicity::NonPeriodic())
Definition: AMReX_FFT_R2X.H:533
R2X(R2X const &)=delete
std::unique_ptr< char, DataDeleter > m_data_1
Definition: AMReX_FFT_R2X.H:102
Plan< T > m_fft_fwd_x
Definition: AMReX_FFT_R2X.H:73
R2X(Box const &domain, Array< std::pair< Boundary, Boundary >, AMREX_SPACEDIM > const &bc, Info const &info=Info{})
Definition: AMReX_FFT_R2X.H:116
std::unique_ptr< MultiBlockCommMetaData > m_cmd_cy2cx
Definition: AMReX_FFT_R2X.H:85
Box m_dom_0
Definition: AMReX_FFT_R2X.H:70
std::unique_ptr< MultiBlockCommMetaData > m_cmd_cx2cy
Definition: AMReX_FFT_R2X.H:80
std::unique_ptr< char, DataDeleter > m_data_2
Definition: AMReX_FFT_R2X.H:103
Plan< T > m_fft_fwd_z
Definition: AMReX_FFT_R2X.H:77
Plan< T > m_fft_bwd_y
Definition: AMReX_FFT_R2X.H:76
cMF m_cx
Definition: AMReX_FFT_R2X.H:98
cMF m_cz
Definition: AMReX_FFT_R2X.H:100
R2X & operator=(R2X const &)=delete
Swap02 m_dtos_z2y
Definition: AMReX_FFT_R2X.H:93
cMF m_cy
Definition: AMReX_FFT_R2X.H:99
Box m_dom_cx
Definition: AMReX_FFT_R2X.H:108
Plan< T > m_fft_bwd_x
Definition: AMReX_FFT_R2X.H:74
Box m_dom_rz
Definition: AMReX_FFT_R2X.H:107
int size() const noexcept
Return the number of FABs in the FabArray.
Definition: AMReX_FabArrayBase.H:109
bool empty() const noexcept
Definition: AMReX_FabArrayBase.H:88
const DistributionMapping & DistributionMap() const noexcept
Return constant reference to associated DistributionMapping.
Definition: AMReX_FabArrayBase.H:130
Box box(int K) const noexcept
Return the Kth Box in the BoxArray. That is, the valid region of the Kth grid.
Definition: AMReX_FabArrayBase.H:100
An Array of FortranArrayBox(FAB)-like Objects.
Definition: AMReX_FabArray.H:344
void define(const BoxArray &bxs, const DistributionMapping &dm, int nvar, int ngrow, const MFInfo &info=MFInfo(), const FabFactory< FAB > &factory=DefaultFabFactory< FAB >())
Define this FabArray identically to that performed by the constructor having an analogous function si...
Definition: AMReX_FabArray.H:2027
A collection (stored as an array) of FArrayBox objects.
Definition: AMReX_MultiFab.H:38
This provides length of period for periodic domains. 0 means it is not periodic in that direction....
Definition: AMReX_Periodicity.H:17
static const Periodicity & NonPeriodic() noexcept
Definition: AMReX_Periodicity.cpp:52
@ FAB
Definition: AMReX_AmrvisConstants.H:86
FA::FABType::value_type * get_fab(FA &fa)
Definition: AMReX_FFT_Helper.H:1303
std::unique_ptr< char, DataDeleter > make_mfs_share(FA1 &fa1, FA2 &fa2)
Definition: AMReX_FFT_Helper.H:1314
DistributionMapping make_iota_distromap(Long n)
Definition: AMReX_FFT.cpp:88
Definition: AMReX_FFT.cpp:7
int MyProcSub() noexcept
my sub-rank in current frame
Definition: AMReX_ParallelContext.H:76
int NProcsSub() noexcept
number of ranks in current frame
Definition: AMReX_ParallelContext.H:74
static int f(amrex::Real t, N_Vector y_data, N_Vector y_rhs, void *user_data)
Definition: AMReX_SundialsIntegrator.H:44
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void swap(T &a, T &b) noexcept
Definition: AMReX_algoim_K.H:113
@ min
Definition: AMReX_ParallelReduce.H:18
std::enable_if_t< std::is_integral_v< T > > ParallelFor(TypeList< CTOs... > ctos, std::array< int, sizeof...(CTOs)> const &runtime_options, T N, F &&f)
Definition: AMReX_CTOParallelForImpl.H:200
BoxND< AMREX_SPACEDIM > Box
Definition: AMReX_BaseFwd.H:27
double second() noexcept
Definition: AMReX_Utility.cpp:922
IntVectND< AMREX_SPACEDIM > IntVect
Definition: AMReX_BaseFwd.H:30
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void ignore_unused(const Ts &...)
This shuts up the compiler about unused variables.
Definition: AMReX.H:111
BoxArray decompose(Box const &domain, int nboxes, Array< bool, AMREX_SPACEDIM > const &decomp={AMREX_D_DECL(true, true, true)}, bool no_overlap=false)
Decompose domain box into BoxArray.
void Abort(const std::string &msg)
Print out message to cerr and exit via abort().
Definition: AMReX.cpp:225
void ParallelCopy(MF &dst, MF const &src, int scomp, int dcomp, int ncomp, IntVect const &ng_src=IntVect(0), IntVect const &ng_dst=IntVect(0), Periodicity const &period=Periodicity::NonPeriodic())
dst = src w/ MPI communication
Definition: AMReX_FabArrayUtility.H:1672
std::array< T, N > Array
Definition: AMReX_Array.H:24
Definition: AMReX_FFT_Helper.H:56
bool batch_mode
Definition: AMReX_FFT_Helper.H:60
int nprocs
Max number of processes to use.
Definition: AMReX_FFT_Helper.H:63
Definition: AMReX_FFT_Helper.H:111
std::conditional_t< std::is_same_v< float, T >, cuComplex, cuDoubleComplex > VendorComplex
Definition: AMReX_FFT_Helper.H:115
Definition: AMReX_FFT_Helper.H:1355
Definition: AMReX_FFT_Helper.H:1378
FabArray memory allocation information.
Definition: AMReX_FabArray.H:66