Block-Structured AMR Software Framework
AMReX_MLPoisson.H
Go to the documentation of this file.
1 #ifndef AMREX_MLPOISSON_H_
2 #define AMREX_MLPOISSON_H_
3 #include <AMReX_Config.H>
4 
5 #include <AMReX_MLCellABecLap.H>
6 #include <AMReX_MLPoisson_K.H>
7 #include <AMReX_MLALaplacian.H>
8 
9 namespace amrex {
10 
11 // del dot grad phi
12 
13 template <typename MF>
15  : public MLCellABecLapT<MF>
16 {
17 public:
18 
19  using FAB = typename MF::fab_type;
20  using RT = typename MF::value_type;
21 
22  using BCType = LinOpBCType;
23  using Location = typename MLLinOpT<MF>::Location;
24 
25  MLPoissonT () = default;
26  MLPoissonT (const Vector<Geometry>& a_geom,
27  const Vector<BoxArray>& a_grids,
28  const Vector<DistributionMapping>& a_dmap,
29  const LPInfo& a_info = LPInfo(),
30  const Vector<FabFactory<FAB> const*>& a_factory = {});
31  MLPoissonT (const Vector<Geometry>& a_geom,
32  const Vector<BoxArray>& a_grids,
33  const Vector<DistributionMapping>& a_dmap,
34  const Vector<iMultiFab const*>& a_overset_mask, // 1: unknown, 0: known
35  const LPInfo& a_info = LPInfo(),
36  const Vector<FabFactory<FAB> const*>& a_factory = {});
37  ~MLPoissonT () override;
38 
39  MLPoissonT (const MLPoissonT<MF>&) = delete;
40  MLPoissonT (MLPoissonT<MF>&&) = delete;
43 
44  void define (const Vector<Geometry>& a_geom,
45  const Vector<BoxArray>& a_grids,
46  const Vector<DistributionMapping>& a_dmap,
47  const LPInfo& a_info = LPInfo(),
48  const Vector<FabFactory<FAB> const*>& a_factory = {});
49 
50  void define (const Vector<Geometry>& a_geom,
51  const Vector<BoxArray>& a_grids,
52  const Vector<DistributionMapping>& a_dmap,
53  const Vector<iMultiFab const*>& a_overset_mask,
54  const LPInfo& a_info = LPInfo(),
55  const Vector<FabFactory<FAB> const*>& a_factory = {});
56 
57  void prepareForSolve () final;
58  [[nodiscard]] bool isSingular (int amrlev) const final { return m_is_singular[amrlev]; }
59  [[nodiscard]] bool isBottomSingular () const final { return m_is_singular[0]; }
60  void Fapply (int amrlev, int mglev, MF& out, const MF& in) const final;
61  void Fsmooth (int amrlev, int mglev, MF& sol, const MF& rhs, int redblack) const final;
62  void FFlux (int amrlev, const MFIter& mfi,
63  const Array<FAB*,AMREX_SPACEDIM>& flux,
64  const FAB& sol, Location loc, int face_only=0) const final;
65 
66  void normalize (int amrlev, int mglev, MF& mf) const final;
67 
68  [[nodiscard]] RT getAScalar () const final { return RT(0.0); }
69  [[nodiscard]] RT getBScalar () const final { return RT(-1.0); }
70  [[nodiscard]] MF const* getACoeffs (int /*amrlev*/, int /*mglev*/) const final { return nullptr; }
71  [[nodiscard]] Array<MF const*,AMREX_SPACEDIM> getBCoeffs (int /*amrlev*/, int /*mglev*/) const final
72  { return {{ AMREX_D_DECL(nullptr,nullptr,nullptr)}}; }
73 
74  [[nodiscard]] std::unique_ptr<MLLinOpT<MF>> makeNLinOp (int grid_size) const final;
75 
76  [[nodiscard]] bool supportNSolve () const final;
77 
78  void copyNSolveSolution (MF& dst, MF const& src) const final;
79 
81  void get_dpdn_on_domain_faces (Array<MF*,AMREX_SPACEDIM> const& dpdn,
82  MF const& phi);
83 
84 private:
85 
87 };
88 
89 template <typename MF>
90 MLPoissonT<MF>::MLPoissonT (const Vector<Geometry>& a_geom,
91  const Vector<BoxArray>& a_grids,
92  const Vector<DistributionMapping>& a_dmap,
93  const LPInfo& a_info,
94  const Vector<FabFactory<FAB> const*>& a_factory)
95 {
96  define(a_geom, a_grids, a_dmap, a_info, a_factory);
97 }
98 
99 template <typename MF>
101  const Vector<BoxArray>& a_grids,
102  const Vector<DistributionMapping>& a_dmap,
103  const Vector<iMultiFab const*>& a_overset_mask,
104  const LPInfo& a_info,
105  const Vector<FabFactory<FAB> const*>& a_factory)
106 {
107  define(a_geom, a_grids, a_dmap, a_overset_mask, a_info, a_factory);
108 }
109 
110 template <typename MF>
111 void
113  const Vector<BoxArray>& a_grids,
114  const Vector<DistributionMapping>& a_dmap,
115  const LPInfo& a_info,
116  const Vector<FabFactory<FAB> const*>& a_factory)
117 {
118  BL_PROFILE("MLPoisson::define()");
119  MLCellABecLapT<MF>::define(a_geom, a_grids, a_dmap, a_info, a_factory);
120 }
121 
122 template <typename MF>
123 void
125  const Vector<BoxArray>& a_grids,
126  const Vector<DistributionMapping>& a_dmap,
127  const Vector<iMultiFab const*>& a_overset_mask,
128  const LPInfo& a_info,
129  const Vector<FabFactory<FAB> const*>& a_factory)
130 {
131  BL_PROFILE("MLPoisson::define(overset)");
132  MLCellABecLapT<MF>::define(a_geom, a_grids, a_dmap, a_overset_mask, a_info, a_factory);
133 }
134 
135 template <typename MF>
136 MLPoissonT<MF>::~MLPoissonT () = default;
137 
138 template <typename MF>
139 void
141 {
142  BL_PROFILE("MLPoisson::prepareForSolve()");
143 
145 
146  m_is_singular.clear();
147  m_is_singular.resize(this->m_num_amr_levels, false);
148  auto itlo = std::find(this->m_lobc[0].begin(), this->m_lobc[0].end(), BCType::Dirichlet);
149  auto ithi = std::find(this->m_hibc[0].begin(), this->m_hibc[0].end(), BCType::Dirichlet);
150  if (itlo == this->m_lobc[0].end() && ithi == this->m_hibc[0].end())
151  { // No Dirichlet
152  for (int alev = 0; alev < this->m_num_amr_levels; ++alev)
153  {
154  // For now this assumes that overset regions are treated as Dirichlet bc's
155  if (this->m_domain_covered[alev] && !this->m_overset_mask[alev][0])
156  {
157  m_is_singular[alev] = true;
158  }
159  }
160  }
161  if (!m_is_singular[0] && this->m_needs_coarse_data_for_bc &&
162  this->m_coarse_fine_bc_type == LinOpBCType::Neumann)
163  {
164  AMREX_ASSERT(this->m_overset_mask[0][0] == nullptr);
165  auto bbox = this->m_grids[0][0].minimalBox();
166  for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
167  if (this->m_lobc[0][idim] == LinOpBCType::Dirichlet) {
168  bbox.growLo(idim,1);
169  }
170  if (this->m_hibc[0][idim] == LinOpBCType::Dirichlet) {
171  bbox.growHi(idim,1);
172  }
173  }
174  if (this->m_geom[0][0].Domain().contains(bbox)) {
175  m_is_singular[0] = true;
176  }
177  }
178 }
179 
180 template <typename MF>
181 void
182 MLPoissonT<MF>::Fapply (int amrlev, int mglev, MF& out, const MF& in) const
183 {
184  BL_PROFILE("MLPoisson::Fapply()");
185 
186  const Real* dxinv = this->m_geom[amrlev][mglev].InvCellSize();
187 
188  AMREX_D_TERM(const RT dhx = RT(dxinv[0]*dxinv[0]);,
189  const RT dhy = RT(dxinv[1]*dxinv[1]);,
190  const RT dhz = RT(dxinv[2]*dxinv[2]););
191 
192 #if (AMREX_SPACEDIM == 3)
193  RT dh0 = this->get_d0(dhx, dhy, dhz);
194  RT dh1 = this->get_d1(dhx, dhy, dhz);
195 #endif
196 
197 #if (AMREX_SPACEDIM < 3)
198  const RT dx = RT(this->m_geom[amrlev][mglev].CellSize(0));
199  const RT probxlo = RT(this->m_geom[amrlev][mglev].ProbLo(0));
200 #endif
201 
202 #ifdef AMREX_USE_GPU
203  if (Gpu::inLaunchRegion() && out.isFusingCandidate() && !this->hasHiddenDimension()) {
204  auto const& xma = in.const_arrays();
205  auto const& yma = out.arrays();
206  if (this->m_overset_mask[amrlev][mglev]) {
208  const auto& osmma = this->m_overset_mask[amrlev][mglev]->const_arrays();
209  ParallelFor(out,
210  [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept
211  {
213  mlpoisson_adotx_os(AMREX_D_DECL(i,j,k), yma[box_no], xma[box_no], osmma[box_no],
214  AMREX_D_DECL(dhx,dhy,dhz));
215  });
216  } else {
217 #if (AMREX_SPACEDIM < 3)
218  if (this->m_has_metric_term) {
219  ParallelFor(out,
220  [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept
221  {
223  mlpoisson_adotx_m(AMREX_D_DECL(i,j,k), yma[box_no], xma[box_no],
224  AMREX_D_DECL(dhx,dhy,dhz), dx, probxlo);
225  });
226  } else
227 #endif
228  {
229  ParallelFor(out,
230  [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept
231  {
233  mlpoisson_adotx(AMREX_D_DECL(i,j,k), yma[box_no], xma[box_no],
234  AMREX_D_DECL(dhx,dhy,dhz));
235  });
236  }
237  }
239  } else
240 #endif
241  {
242 #ifdef AMREX_USE_OMP
243 #pragma omp parallel if (Gpu::notInLaunchRegion())
244 #endif
245  for (MFIter mfi(out, TilingIfNotGPU()); mfi.isValid(); ++mfi)
246  {
247  const Box& bx = mfi.tilebox();
248  const auto& xfab = in.array(mfi);
249  const auto& yfab = out.array(mfi);
250 
251  if (this->m_overset_mask[amrlev][mglev]) {
253  const auto& osm = this->m_overset_mask[amrlev][mglev]->const_array(mfi);
255  {
257  mlpoisson_adotx_os(AMREX_D_DECL(i,j,k), yfab, xfab, osm,
258  AMREX_D_DECL(dhx,dhy,dhz));
259  });
260  } else {
261 #if (AMREX_SPACEDIM == 3)
262  if (this->hasHiddenDimension()) {
263  Box const& bx2d = this->compactify(bx);
264  const auto& xfab2d = this->compactify(xfab);
265  const auto& yfab2d = this->compactify(yfab);
266  AMREX_HOST_DEVICE_PARALLEL_FOR_3D(bx2d, i, j, k,
267  {
269  TwoD::mlpoisson_adotx(i, j, yfab2d, xfab2d, dh0, dh1);
270  });
271  } else {
273  {
274  mlpoisson_adotx(i, j, k, yfab, xfab, dhx, dhy, dhz);
275  });
276  }
277 #elif (AMREX_SPACEDIM == 2)
278  if (this->m_has_metric_term) {
280  {
282  mlpoisson_adotx_m(i, j, yfab, xfab, dhx, dhy, dx, probxlo);
283  });
284  } else {
286  {
288  mlpoisson_adotx(i, j, yfab, xfab, dhx, dhy);
289  });
290  }
291 #elif (AMREX_SPACEDIM == 1)
292  if (this->m_has_metric_term) {
294  {
296  mlpoisson_adotx_m(i, yfab, xfab, dhx, dx, probxlo);
297  });
298  } else {
300  {
302  mlpoisson_adotx(i, yfab, xfab, dhx);
303  });
304  }
305 #endif
306  }
307  }
308  }
309 }
310 
311 template <typename MF>
312 void
313 MLPoissonT<MF>::normalize (int amrlev, int mglev, MF& mf) const
314 {
315  amrex::ignore_unused(amrlev,mglev,mf);
316 #if (AMREX_SPACEDIM != 3)
317  BL_PROFILE("MLPoisson::normalize()");
318 
319  if (!this->m_has_metric_term) { return; }
320 
321  const Real* dxinv = this->m_geom[amrlev][mglev].InvCellSize();
322  AMREX_D_TERM(const RT dhx = RT(dxinv[0]*dxinv[0]);,
323  const RT dhy = RT(dxinv[1]*dxinv[1]);,
324  const RT dhz = RT(dxinv[2]*dxinv[2]););
325  const RT dx = RT(this->m_geom[amrlev][mglev].CellSize(0));
326  const RT probxlo = RT(this->m_geom[amrlev][mglev].ProbLo(0));
327 
328 #ifdef AMREX_USE_GPU
329  if (Gpu::inLaunchRegion() && mf.isFusingCandidate()) {
330  auto const& ma = mf.arrays();
331  ParallelFor(mf,
332  [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept
333  {
334  mlpoisson_normalize(i,j,k, ma[box_no], AMREX_D_DECL(dhx,dhy,dhz), dx, probxlo);
335  });
337  } else
338 #endif
339  {
340 #ifdef AMREX_USE_OMP
341 #pragma omp parallel if (Gpu::notInLaunchRegion())
342 #endif
343  for (MFIter mfi(mf, TilingIfNotGPU()); mfi.isValid(); ++mfi)
344  {
345  const Box& bx = mfi.tilebox();
346  const auto& fab = mf.array(mfi);
347 
348 #if (AMREX_SPACEDIM == 2)
350  {
351  mlpoisson_normalize(i,j,k, fab, dhx, dhy, dx, probxlo);
352  });
353 #else
355  {
356  mlpoisson_normalize(i,j,k, fab, dhx, dx, probxlo);
357  });
358 #endif
359  }
360  }
361 #endif
362 }
363 
364 template <typename MF>
365 void
366 MLPoissonT<MF>::Fsmooth (int amrlev, int mglev, MF& sol, const MF& rhs, int redblack) const
367 {
368  BL_PROFILE("MLPoisson::Fsmooth()");
369 
370  MF Ax;
371  if (! this->m_use_gauss_seidel) { // jacobi
372  Ax.define(sol.boxArray(), sol.DistributionMap(), sol.nComp(), 0);
373  Fapply(amrlev, mglev, Ax, sol);
374  }
375 
376  const auto& undrrelxr = this->m_undrrelxr[amrlev][mglev];
377  const auto& maskvals = this->m_maskvals [amrlev][mglev];
378 
379  OrientationIter oitr;
380 
381  const auto& f0 = undrrelxr[oitr()]; ++oitr;
382  const auto& f1 = undrrelxr[oitr()]; ++oitr;
383 #if (AMREX_SPACEDIM > 1)
384  const auto& f2 = undrrelxr[oitr()]; ++oitr;
385  const auto& f3 = undrrelxr[oitr()]; ++oitr;
386 #if (AMREX_SPACEDIM > 2)
387  const auto& f4 = undrrelxr[oitr()]; ++oitr;
388  const auto& f5 = undrrelxr[oitr()]; ++oitr;
389 #endif
390 #endif
391 
392  const MultiMask& mm0 = maskvals[0];
393  const MultiMask& mm1 = maskvals[1];
394 #if (AMREX_SPACEDIM > 1)
395  const MultiMask& mm2 = maskvals[2];
396  const MultiMask& mm3 = maskvals[3];
397 #if (AMREX_SPACEDIM > 2)
398  const MultiMask& mm4 = maskvals[4];
399  const MultiMask& mm5 = maskvals[5];
400 #endif
401 #endif
402 
403  const Real* dxinv = this->m_geom[amrlev][mglev].InvCellSize();
404  AMREX_D_TERM(const RT dhx = RT(dxinv[0]*dxinv[0]);,
405  const RT dhy = RT(dxinv[1]*dxinv[1]);,
406  const RT dhz = RT(dxinv[2]*dxinv[2]););
407 
408 #if (AMREX_SPACEDIM == 3)
409  RT dh0 = RT(this->get_d0(dhx, dhy, dhz));
410  RT dh1 = RT(this->get_d1(dhx, dhy, dhz));
411 #endif
412 
413 #if (AMREX_SPACEDIM < 3)
414  const RT dx = RT(this->m_geom[amrlev][mglev].CellSize(0));
415  const RT probxlo = RT(this->m_geom[amrlev][mglev].ProbLo(0));
416 #endif
417 
418  MFItInfo mfi_info;
419  if (Gpu::notInLaunchRegion()) { mfi_info.EnableTiling().SetDynamic(true); }
420 
421 #ifdef AMREX_USE_GPU
422  if (Gpu::inLaunchRegion() && sol.isFusingCandidate()
423  && ! this->hasHiddenDimension())
424  {
425  const auto& m0ma = mm0.const_arrays();
426  const auto& m1ma = mm1.const_arrays();
427 #if (AMREX_SPACEDIM > 1)
428  const auto& m2ma = mm2.const_arrays();
429  const auto& m3ma = mm3.const_arrays();
430 #if (AMREX_SPACEDIM > 2)
431  const auto& m4ma = mm4.const_arrays();
432  const auto& m5ma = mm5.const_arrays();
433 #endif
434 #endif
435 
436  const auto& solnma = sol.arrays();
437  const auto& rhsma = rhs.const_arrays();
438 
439  AMREX_ALWAYS_ASSERT(rhs.nGrowVect() == 0);
440 
441  const auto& f0ma = f0.const_arrays();
442  const auto& f1ma = f1.const_arrays();
443 #if (AMREX_SPACEDIM > 1)
444  const auto& f2ma = f2.const_arrays();
445  const auto& f3ma = f3.const_arrays();
446 #if (AMREX_SPACEDIM > 2)
447  const auto& f4ma = f4.const_arrays();
448  const auto& f5ma = f5.const_arrays();
449 #endif
450 #endif
451 
452  if (this->m_overset_mask[amrlev][mglev]) {
454  const auto& osmma = this->m_overset_mask[amrlev][mglev]->const_arrays();
455  if (this->m_use_gauss_seidel) {
456  ParallelFor(sol,
457  [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept
458  {
459  Box vbx(rhsma[box_no]);
460  mlpoisson_gsrb_os(i, j, k, solnma[box_no], rhsma[box_no],
461  osmma[box_no], AMREX_D_DECL(dhx, dhy, dhz),
462  f0ma[box_no], m0ma[box_no],
463  f1ma[box_no], m1ma[box_no],
464 #if (AMREX_SPACEDIM > 1)
465  f2ma[box_no], m2ma[box_no],
466  f3ma[box_no], m3ma[box_no],
467 #if (AMREX_SPACEDIM > 2)
468  f4ma[box_no], m4ma[box_no],
469  f5ma[box_no], m5ma[box_no],
470 #endif
471 #endif
472  vbx, redblack);
473  });
474  } else {
475  const auto& axma = Ax.const_arrays();
476  ParallelFor(sol,
477  [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept
478  {
479  Box vbx(rhsma[box_no]);
480  mlpoisson_jacobi_os(i, j, k, solnma[box_no], rhsma[box_no],
481  axma[box_no], osmma[box_no],
482  AMREX_D_DECL(dhx, dhy, dhz),
483  f0ma[box_no], m0ma[box_no],
484  f1ma[box_no], m1ma[box_no],
485 #if (AMREX_SPACEDIM > 1)
486  f2ma[box_no], m2ma[box_no],
487  f3ma[box_no], m3ma[box_no],
488 #if (AMREX_SPACEDIM > 2)
489  f4ma[box_no], m4ma[box_no],
490  f5ma[box_no], m5ma[box_no],
491 #endif
492 #endif
493  vbx);
494  });
495  }
496  }
497 #if (AMREX_SPACEDIM < 3)
498  else if (this->m_has_metric_term) {
499  if (this->m_use_gauss_seidel) {
500  ParallelFor(sol,
501  [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept
502  {
503  Box vbx(rhsma[box_no]);
504  mlpoisson_gsrb_m(i, j, k, solnma[box_no], rhsma[box_no],
505  AMREX_D_DECL(dhx, dhy, dhz),
506  f0ma[box_no], m0ma[box_no],
507  f1ma[box_no], m1ma[box_no],
508 #if (AMREX_SPACEDIM > 1)
509  f2ma[box_no], m2ma[box_no],
510  f3ma[box_no], m3ma[box_no],
511 #endif
512  vbx, redblack,
513  dx, probxlo);
514  });
515  } else {
516  const auto& axma = Ax.const_arrays();
517  ParallelFor(sol,
518  [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept
519  {
520  Box vbx(rhsma[box_no]);
521  mlpoisson_jacobi_m(i, j, k, solnma[box_no], rhsma[box_no],
522  axma[box_no], AMREX_D_DECL(dhx, dhy, dhz),
523  f0ma[box_no], m0ma[box_no],
524  f1ma[box_no], m1ma[box_no],
525 #if (AMREX_SPACEDIM > 1)
526  f2ma[box_no], m2ma[box_no],
527  f3ma[box_no], m3ma[box_no],
528 #endif
529  vbx, dx, probxlo);
530  });
531  }
532  }
533 #endif
534  else {
535  if (this->m_use_gauss_seidel) {
536  ParallelFor(sol,
537  [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept
538  {
539  Box vbx(rhsma[box_no]);
540  mlpoisson_gsrb(i, j, k, solnma[box_no], rhsma[box_no],
541  AMREX_D_DECL(dhx, dhy, dhz),
542  f0ma[box_no], m0ma[box_no],
543  f1ma[box_no], m1ma[box_no],
544 #if (AMREX_SPACEDIM > 1)
545  f2ma[box_no], m2ma[box_no],
546  f3ma[box_no], m3ma[box_no],
547 #if (AMREX_SPACEDIM > 2)
548  f4ma[box_no], m4ma[box_no],
549  f5ma[box_no], m5ma[box_no],
550 #endif
551 #endif
552  vbx, redblack);
553  });
554  } else {
555  const auto& axma = Ax.const_arrays();
556  ParallelFor(sol,
557  [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept
558  {
559  Box vbx(rhsma[box_no]);
560  mlpoisson_jacobi(i, j, k, solnma[box_no], rhsma[box_no],
561  axma[box_no], AMREX_D_DECL(dhx, dhy, dhz),
562  f0ma[box_no], m0ma[box_no],
563  f1ma[box_no], m1ma[box_no],
564 #if (AMREX_SPACEDIM > 1)
565  f2ma[box_no], m2ma[box_no],
566  f3ma[box_no], m3ma[box_no],
567 #if (AMREX_SPACEDIM > 2)
568  f4ma[box_no], m4ma[box_no],
569  f5ma[box_no], m5ma[box_no],
570 #endif
571 #endif
572  vbx);
573  });
574  }
575  }
576  } else
577 #endif
578  {
579 #ifdef AMREX_USE_OMP
580 #pragma omp parallel if (Gpu::notInLaunchRegion())
581 #endif
582  for (MFIter mfi(sol,mfi_info); mfi.isValid(); ++mfi)
583  {
584  const auto& m0 = mm0.array(mfi);
585  const auto& m1 = mm1.array(mfi);
586 #if (AMREX_SPACEDIM > 1)
587  const auto& m2 = mm2.array(mfi);
588  const auto& m3 = mm3.array(mfi);
589 #if (AMREX_SPACEDIM > 2)
590  const auto& m4 = mm4.array(mfi);
591  const auto& m5 = mm5.array(mfi);
592 #endif
593 #endif
594 
595  const Box& tbx = mfi.tilebox();
596  const Box& vbx = mfi.validbox();
597  const auto& solnfab = sol.array(mfi);
598  const auto& rhsfab = rhs.array(mfi);
599 
600  const auto& f0fab = f0.array(mfi);
601  const auto& f1fab = f1.array(mfi);
602 #if (AMREX_SPACEDIM > 1)
603  const auto& f2fab = f2.array(mfi);
604  const auto& f3fab = f3.array(mfi);
605 #if (AMREX_SPACEDIM > 2)
606  const auto& f4fab = f4.array(mfi);
607  const auto& f5fab = f5.array(mfi);
608 #endif
609 #endif
610 
611 #if (AMREX_SPACEDIM == 1)
612  if (this->m_overset_mask[amrlev][mglev]) {
614  const auto& osm = this->m_overset_mask[amrlev][mglev]->const_array(mfi);
615  if (this->m_use_gauss_seidel) {
616  AMREX_HOST_DEVICE_PARALLEL_FOR_3D ( tbx, i, j, k,
617  {
618  mlpoisson_gsrb_os(i, j, k, solnfab, rhsfab, osm, dhx,
619  f0fab, m0,
620  f1fab, m1,
621  vbx, redblack);
622  });
623  } else {
624  const auto& axfab = Ax.const_array(mfi);
625  AMREX_HOST_DEVICE_PARALLEL_FOR_3D ( tbx, i, j, k,
626  {
627  mlpoisson_jacobi_os(i, j, k, solnfab, rhsfab, axfab,
628  osm, dhx,
629  f0fab, m0,
630  f1fab, m1,
631  vbx);
632  });
633  }
634  } else if (this->m_has_metric_term) {
635  if (this->m_use_gauss_seidel) {
636  AMREX_HOST_DEVICE_PARALLEL_FOR_3D ( tbx, i, j, k,
637  {
638  mlpoisson_gsrb_m(i, j, k, solnfab, rhsfab, dhx,
639  f0fab, m0,
640  f1fab, m1,
641  vbx, redblack,
642  dx, probxlo);
643  });
644  } else {
645  const auto& axfab = Ax.const_array(mfi);
646  AMREX_HOST_DEVICE_PARALLEL_FOR_3D ( tbx, i, j, k,
647  {
648  mlpoisson_jacobi_m(i, j, k, solnfab, rhsfab, axfab, dhx,
649  f0fab, m0,
650  f1fab, m1,
651  vbx, dx, probxlo);
652  });
653  }
654  } else {
655  if (this->m_use_gauss_seidel) {
656  AMREX_HOST_DEVICE_PARALLEL_FOR_3D ( tbx, i, j, k,
657  {
658  mlpoisson_gsrb(i, j, k, solnfab, rhsfab, dhx,
659  f0fab, m0,
660  f1fab, m1,
661  vbx, redblack);
662  });
663  } else {
664  const auto& axfab = Ax.const_array(mfi);
665  AMREX_HOST_DEVICE_PARALLEL_FOR_3D ( tbx, i, j, k,
666  {
667  mlpoisson_jacobi(i, j, k, solnfab, rhsfab, axfab, dhx,
668  f0fab, m0,
669  f1fab, m1,
670  vbx);
671  });
672  }
673  }
674 #endif
675 
676 #if (AMREX_SPACEDIM == 2)
677  if (this->m_overset_mask[amrlev][mglev]) {
679  const auto& osm = this->m_overset_mask[amrlev][mglev]->const_array(mfi);
680  if (this->m_use_gauss_seidel) {
681  AMREX_HOST_DEVICE_PARALLEL_FOR_3D ( tbx, i, j, k,
682  {
683  mlpoisson_gsrb_os(i, j, k, solnfab, rhsfab, osm, dhx, dhy,
684  f0fab, m0,
685  f1fab, m1,
686  f2fab, m2,
687  f3fab, m3,
688  vbx, redblack);
689  });
690  } else {
691  const auto& axfab = Ax.const_array(mfi);
692  AMREX_HOST_DEVICE_PARALLEL_FOR_3D ( tbx, i, j, k,
693  {
694  mlpoisson_jacobi_os(i, j, k, solnfab, rhsfab, axfab,
695  osm, dhx, dhy,
696  f0fab, m0,
697  f1fab, m1,
698  f2fab, m2,
699  f3fab, m3,
700  vbx);
701  });
702  }
703  } else if (this->m_has_metric_term) {
704  if (this->m_use_gauss_seidel) {
705  AMREX_HOST_DEVICE_PARALLEL_FOR_3D ( tbx, i, j, k,
706  {
707  mlpoisson_gsrb_m(i, j, k, solnfab, rhsfab, dhx, dhy,
708  f0fab, m0,
709  f1fab, m1,
710  f2fab, m2,
711  f3fab, m3,
712  vbx, redblack,
713  dx, probxlo);
714  });
715  } else {
716  const auto& axfab = Ax.const_array(mfi);
717  AMREX_HOST_DEVICE_PARALLEL_FOR_3D ( tbx, i, j, k,
718  {
719  mlpoisson_jacobi_m(i, j, k, solnfab, rhsfab, axfab, dhx, dhy,
720  f0fab, m0,
721  f1fab, m1,
722  f2fab, m2,
723  f3fab, m3,
724  vbx, dx, probxlo);
725  });
726  }
727  } else {
728  if (this->m_use_gauss_seidel) {
729  AMREX_HOST_DEVICE_PARALLEL_FOR_3D ( tbx, i, j, k,
730  {
731  mlpoisson_gsrb(i, j, k, solnfab, rhsfab, dhx, dhy,
732  f0fab, m0,
733  f1fab, m1,
734  f2fab, m2,
735  f3fab, m3,
736  vbx, redblack);
737  });
738  } else {
739  const auto& axfab = Ax.const_array(mfi);
740  AMREX_HOST_DEVICE_PARALLEL_FOR_3D ( tbx, i, j, k,
741  {
742  mlpoisson_jacobi(i, j, k, solnfab, rhsfab, axfab, dhx, dhy,
743  f0fab, m0,
744  f1fab, m1,
745  f2fab, m2,
746  f3fab, m3,
747  vbx);
748  });
749  }
750  }
751 #endif
752 
753 #if (AMREX_SPACEDIM == 3)
754  if (this->m_overset_mask[amrlev][mglev]) {
756  const auto& osm = this->m_overset_mask[amrlev][mglev]->const_array(mfi);
757  if (this->m_use_gauss_seidel) {
758  AMREX_HOST_DEVICE_PARALLEL_FOR_3D ( tbx, i, j, k,
759  {
760  mlpoisson_gsrb_os(i, j, k, solnfab, rhsfab, osm, dhx, dhy, dhz,
761  f0fab, m0,
762  f1fab, m1,
763  f2fab, m2,
764  f3fab, m3,
765  f4fab, m4,
766  f5fab, m5,
767  vbx, redblack);
768  });
769  } else {
770  const auto& axfab = Ax.const_array(mfi);
771  AMREX_HOST_DEVICE_PARALLEL_FOR_3D ( tbx, i, j, k,
772  {
773  mlpoisson_jacobi_os(i, j, k, solnfab, rhsfab, axfab,
774  osm, dhx, dhy, dhz,
775  f0fab, m0,
776  f1fab, m1,
777  f2fab, m2,
778  f3fab, m3,
779  f4fab, m4,
780  f5fab, m5,
781  vbx);
782  });
783  }
784  } else if (this->hasHiddenDimension()) {
785  Box const& tbx_2d = this->compactify(tbx);
786  Box const& vbx_2d = this->compactify(vbx);
787  const auto& solnfab_2d = this->compactify(solnfab);
788  const auto& rhsfab_2d = this->compactify(rhsfab);
789  const auto& f0fab_2d = this->compactify(this->get_d0(f0fab,f1fab,f2fab));
790  const auto& f1fab_2d = this->compactify(this->get_d1(f0fab,f1fab,f2fab));
791  const auto& f2fab_2d = this->compactify(this->get_d0(f3fab,f4fab,f5fab));
792  const auto& f3fab_2d = this->compactify(this->get_d1(f3fab,f4fab,f5fab));
793  const auto& m0_2d = this->compactify(this->get_d0(m0,m1,m2));
794  const auto& m1_2d = this->compactify(this->get_d1(m0,m1,m2));
795  const auto& m2_2d = this->compactify(this->get_d0(m3,m4,m5));
796  const auto& m3_2d = this->compactify(this->get_d1(m3,m4,m5));
797  if (this->m_use_gauss_seidel) {
798  AMREX_HOST_DEVICE_PARALLEL_FOR_3D ( tbx_2d, i, j, k,
799  {
800  TwoD::mlpoisson_gsrb(i, j, k, solnfab_2d, rhsfab_2d, dh0, dh1,
801  f0fab_2d, m0_2d,
802  f1fab_2d, m1_2d,
803  f2fab_2d, m2_2d,
804  f3fab_2d, m3_2d,
805  vbx_2d, redblack);
806  });
807  } else {
808  const auto& axfab = Ax.const_array(mfi);
809  const auto& axfab_2d = this->compactify(axfab);
810  AMREX_HOST_DEVICE_PARALLEL_FOR_3D ( tbx_2d, i, j, k,
811  {
812  TwoD::mlpoisson_jacobi(i, j, k, solnfab_2d, rhsfab_2d,
813  axfab_2d, dh0, dh1,
814  f0fab_2d, m0_2d,
815  f1fab_2d, m1_2d,
816  f2fab_2d, m2_2d,
817  f3fab_2d, m3_2d,
818  vbx_2d);
819  });
820  }
821  } else {
822  if (this->m_use_gauss_seidel) {
823  AMREX_HOST_DEVICE_PARALLEL_FOR_3D ( tbx, i, j, k,
824  {
825  mlpoisson_gsrb(i, j, k, solnfab, rhsfab, dhx, dhy, dhz,
826  f0fab, m0,
827  f1fab, m1,
828  f2fab, m2,
829  f3fab, m3,
830  f4fab, m4,
831  f5fab, m5,
832  vbx, redblack);
833  });
834  } else {
835  const auto& axfab = Ax.const_array(mfi);
836  AMREX_HOST_DEVICE_PARALLEL_FOR_3D ( tbx, i, j, k,
837  {
838  mlpoisson_jacobi(i, j, k, solnfab, rhsfab, axfab,
839  dhx, dhy, dhz,
840  f0fab, m0,
841  f1fab, m1,
842  f2fab, m2,
843  f3fab, m3,
844  f4fab, m4,
845  f5fab, m5,
846  vbx);
847  });
848  }
849  }
850 #endif
851  }
852  }
853 }
854 
855 template <typename MF>
856 void
857 MLPoissonT<MF>::FFlux (int amrlev, const MFIter& mfi,
858  const Array<FAB*,AMREX_SPACEDIM>& flux,
859  const FAB& sol, Location, const int face_only) const
860 {
861  BL_PROFILE("MLPoisson::FFlux()");
862 
863  const int mglev = 0;
864  const Box& box = mfi.tilebox();
865  const Real* dxinv = this->m_geom[amrlev][mglev].InvCellSize();
866 
867  AMREX_D_TERM(const auto& fxarr = flux[0]->array();,
868  const auto& fyarr = flux[1]->array();,
869  const auto& fzarr = flux[2]->array(););
870  const auto& solarr = sol.array();
871 
872 #if (AMREX_SPACEDIM != 3)
873  const RT dx = RT(this->m_geom[amrlev][mglev].CellSize(0));
874  const RT probxlo = RT(this->m_geom[amrlev][mglev].ProbLo(0));
875 #endif
876 
877 #if (AMREX_SPACEDIM == 3)
878  if (face_only) {
879  if (this->hiddenDirection() != 0) {
880  RT fac = RT(dxinv[0]);
881  Box blo = amrex::bdryLo(box, 0);
882  int blen = box.length(0);
884  {
885  mlpoisson_flux_xface(tbox, fxarr, solarr, fac, blen);
886  });
887  } else {
888  flux[0]->template setVal<RunOn::Device>(RT(0.0));
889  }
890  if (this->hiddenDirection() != 1) {
891  RT fac = RT(dxinv[1]);
892  Box blo = amrex::bdryLo(box, 1);
893  int blen = box.length(1);
895  {
896  mlpoisson_flux_yface(tbox, fyarr, solarr, fac, blen);
897  });
898  } else {
899  flux[1]->template setVal<RunOn::Device>(RT(0.0));
900  }
901  if (this->hiddenDirection() != 2) {
902  RT fac = RT(dxinv[2]);
903  Box blo = amrex::bdryLo(box, 2);
904  int blen = box.length(2);
906  {
907  mlpoisson_flux_zface(tbox, fzarr, solarr, fac, blen);
908  });
909  } else {
910  flux[2]->template setVal<RunOn::Device>(RT(0.0));
911  }
912  } else {
913  if (this->hiddenDirection() != 0) {
914  RT fac = RT(dxinv[0]);
915  Box bflux = amrex::surroundingNodes(box, 0);
916  AMREX_LAUNCH_HOST_DEVICE_LAMBDA ( bflux, tbox,
917  {
918  mlpoisson_flux_x(tbox, fxarr, solarr, fac);
919  });
920  } else {
921  flux[0]->template setVal<RunOn::Device>(RT(0.0));
922  }
923  if (this->hiddenDirection() != 1) {
924  RT fac = RT(dxinv[1]);
925  Box bflux = amrex::surroundingNodes(box, 1);
926  AMREX_LAUNCH_HOST_DEVICE_LAMBDA ( bflux, tbox,
927  {
928  mlpoisson_flux_y(tbox, fyarr, solarr, fac);
929  });
930  } else {
931  flux[1]->template setVal<RunOn::Device>(RT(0.0));
932  }
933  if (this->hiddenDirection() != 2) {
934  RT fac = RT(dxinv[2]);
935  Box bflux = amrex::surroundingNodes(box, 2);
936  AMREX_LAUNCH_HOST_DEVICE_LAMBDA ( bflux, tbox,
937  {
938  mlpoisson_flux_z(tbox, fzarr, solarr, fac);
939  });
940  } else {
941  flux[2]->template setVal<RunOn::Device>(RT(0.0));
942  }
943  }
944 #elif (AMREX_SPACEDIM == 2)
945  if (face_only) {
946  if (this->hiddenDirection() != 0) {
947  RT fac = RT(dxinv[0]);
948  Box blo = amrex::bdryLo(box, 0);
949  int blen = box.length(0);
950  if (this->m_has_metric_term) {
952  {
953  mlpoisson_flux_xface_m(tbox, fxarr, solarr, fac, blen, dx, probxlo);
954  });
955  } else {
957  {
958  mlpoisson_flux_xface(tbox, fxarr, solarr, fac, blen);
959  });
960  }
961  } else {
962  flux[0]->template setVal<RunOn::Device>(RT(0.0));
963  }
964  if (this->hiddenDirection() != 1) {
965  RT fac = RT(dxinv[1]);
966  Box blo = amrex::bdryLo(box, 1);
967  int blen = box.length(1);
968  if (this->m_has_metric_term) {
970  {
971  mlpoisson_flux_yface_m(tbox, fyarr, solarr, fac, blen, dx, probxlo);
972  });
973  } else {
975  {
976  mlpoisson_flux_yface(tbox, fyarr, solarr, fac, blen);
977  });
978  }
979  } else {
980  flux[1]->template setVal<RunOn::Device>(RT(0.0));
981  }
982  } else {
983  if (this->hiddenDirection() != 0) {
984  RT fac = RT(dxinv[0]);
985  Box bflux = amrex::surroundingNodes(box, 0);
986  if (this->m_has_metric_term) {
987  AMREX_LAUNCH_HOST_DEVICE_LAMBDA ( bflux, tbox,
988  {
989  mlpoisson_flux_x_m(tbox, fxarr, solarr, fac, dx, probxlo);
990  });
991  } else {
992  AMREX_LAUNCH_HOST_DEVICE_LAMBDA ( bflux, tbox,
993  {
994  mlpoisson_flux_x(tbox, fxarr, solarr, fac);
995  });
996  }
997  } else {
998  flux[0]->template setVal<RunOn::Device>(RT(0.0));
999  }
1000  if (this->hiddenDirection() != 1) {
1001  RT fac = RT(dxinv[1]);
1002  Box bflux = amrex::surroundingNodes(box, 1);
1003  if (this->m_has_metric_term) {
1004  AMREX_LAUNCH_HOST_DEVICE_LAMBDA ( bflux, tbox,
1005  {
1006  mlpoisson_flux_y_m(tbox, fyarr, solarr, fac, dx, probxlo);
1007  });
1008  } else {
1009  AMREX_LAUNCH_HOST_DEVICE_LAMBDA ( bflux, tbox,
1010  {
1011  mlpoisson_flux_y(tbox, fyarr, solarr, fac);
1012  });
1013  }
1014  } else {
1015  flux[1]->template setVal<RunOn::Device>(RT(0.0));
1016  }
1017  }
1018 #else
1019  if (face_only) {
1020  RT fac = RT(dxinv[0]);
1021  Box blo = amrex::bdryLo(box, 0);
1022  int blen = box.length(0);
1023  if (this->m_has_metric_term) {
1024  AMREX_LAUNCH_HOST_DEVICE_LAMBDA ( blo, tbox,
1025  {
1026  mlpoisson_flux_xface_m(tbox, fxarr, solarr, fac, blen, dx, probxlo);
1027  });
1028  } else {
1029  AMREX_LAUNCH_HOST_DEVICE_LAMBDA ( blo, tbox,
1030  {
1031  mlpoisson_flux_xface(tbox, fxarr, solarr, fac, blen);
1032  });
1033  }
1034  } else {
1035  RT fac = RT(dxinv[0]);
1036  Box bflux = amrex::surroundingNodes(box, 0);
1037  if (this->m_has_metric_term) {
1038  AMREX_LAUNCH_HOST_DEVICE_LAMBDA ( bflux, tbox,
1039  {
1040  mlpoisson_flux_x_m(tbox, fxarr, solarr, fac, dx, probxlo);
1041  });
1042  } else {
1043  AMREX_LAUNCH_HOST_DEVICE_LAMBDA ( bflux, tbox,
1044  {
1045  mlpoisson_flux_x(tbox, fxarr, solarr, fac);
1046  });
1047  }
1048  }
1049 #endif
1050 }
1051 
1052 template <typename MF>
1053 bool
1055 {
1056  bool support = true;
1057  if (this->m_domain_covered[0]) { support = false; }
1058  if (this->doAgglomeration()) { support = false; }
1059  if (AMREX_SPACEDIM != 3) { support = false; }
1060  return support;
1061 }
1062 
1063 template <typename MF>
1064 std::unique_ptr<MLLinOpT<MF>>
1065 MLPoissonT<MF>::makeNLinOp (int grid_size) const
1066 {
1067  const Geometry& geom = this->m_geom[0].back();
1068  const BoxArray& ba = this->makeNGrids(grid_size);
1069 
1071  {
1072  const std::vector<std::vector<int> >& sfc = DistributionMapping::makeSFC(ba);
1073  Vector<int> pmap(ba.size());
1075  const int nprocs = ParallelDescriptor::NProcs();
1076  for (int iproc = 0; iproc < nprocs; ++iproc) {
1077  for (int ibox : sfc[iproc]) {
1078  pmap[ibox] = iproc;
1079  }
1080  }
1081  dm.define(std::move(pmap));
1082  }
1083 
1084  LPInfo minfo{};
1085  minfo.has_metric_term = this->info.has_metric_term;
1086 
1087  std::unique_ptr<MLLinOpT<MF>> r{new MLALaplacianT<MF>({geom}, {ba}, {dm}, minfo)};
1088  auto nop = dynamic_cast<MLALaplacianT<MF>*>(r.get());
1089  if (!nop) {
1090  return nullptr;
1091  }
1092 
1093  nop->m_parent = this;
1094 
1095  nop->setMaxOrder(this->maxorder);
1096  nop->setVerbose(this->verbose);
1097 
1098  nop->setDomainBC(this->m_lobc, this->m_hibc);
1099 
1100  if (this->needsCoarseDataForBC())
1101  {
1102  const Real* dx0 = this->m_geom[0][0].CellSize();
1104  fac *= Real(0.5);
1105  RealVect cbloc {AMREX_D_DECL(dx0[0]*fac[0], dx0[1]*fac[1], dx0[2]*fac[2])};
1106  nop->setCoarseFineBCLocation(cbloc);
1107  }
1108 
1109  nop->setScalars(1.0, -1.0);
1110 
1111  const Real* dxinv = geom.InvCellSize();
1112  RT dxscale = RT(dxinv[0]);
1113 #if (AMREX_SPACEDIM >= 2)
1114  dxscale = std::max(dxscale,RT(dxinv[1]));
1115 #endif
1116 #if (AMREX_SPACEDIM == 3)
1117  dxscale = std::max(dxscale,RT(dxinv[2]));
1118 #endif
1119 
1120  MF alpha(ba, dm, 1, 0);
1121  alpha.setVal(RT(1.e30)*dxscale*dxscale);
1122 
1123  MF foo(this->m_grids[0].back(), this->m_dmap[0].back(), 1, 0, MFInfo().SetAlloc(false));
1124  const FabArrayBase::CPC& cpc = alpha.getCPC(IntVect(0),foo,IntVect(0),Periodicity::NonPeriodic());
1125  alpha.setVal(RT(0.0), cpc, 0, 1);
1126 
1127  nop->setACoeffs(0, alpha);
1128 
1129  return r;
1130 }
1131 
1132 template <typename MF>
1133 void
1134 MLPoissonT<MF>::copyNSolveSolution (MF& dst, MF const& src) const
1135 {
1136  dst.ParallelCopy(src);
1137 }
1138 
1139 template <typename MF>
1140 void
1142  MF const& phi)
1143 {
1144  BL_PROFILE("MLPoisson::dpdn_faces()");
1145 
1146  // We do not need to call applyBC because this function is used by the
1147  // OpenBC solver after solver has converged. That means the BC has been
1148  // filled to check the residual.
1149 
1150  Box const& domain0 = this->m_geom[0][0].Domain();
1151  AMREX_D_TERM(const RT dxi = RT(this->m_geom[0][0].InvCellSize(0));,
1152  const RT dyi = RT(this->m_geom[0][0].InvCellSize(1));,
1153  const RT dzi = RT(this->m_geom[0][0].InvCellSize(2));)
1154 
1155 #ifdef AMREX_USE_OMP
1156 #pragma omp parallel if (Gpu::notInLaunchRegion())
1157 #endif
1158  for (MFIter mfi(phi); mfi.isValid(); ++mfi)
1159  {
1160  Box const& vbx = mfi.validbox();
1161  for (OrientationIter oit; oit.isValid(); ++oit) {
1162  Orientation face = oit();
1163  if (vbx[face] == domain0[face]) {
1164  int dir = face.coordDir();
1165  auto const& p = phi.const_array(mfi);
1166  auto const& gp = dpdn[dir]->array(mfi);
1167  Box const& b2d = amrex::bdryNode(vbx,face);
1168  if (dir == 0) {
1169  // because it's dphi/dn, not dphi/dx.
1170  RT fac = dxi * (face.isLow() ? RT(-1.0) : RT(1.));
1171  AMREX_HOST_DEVICE_PARALLEL_FOR_3D(b2d, i, j, k,
1172  {
1173  gp(i,j,k) = fac * (p(i,j,k) - p(i-1,j,k));
1174  });
1175  }
1176 #if (AMREX_SPACEDIM > 1)
1177  else if (dir == 1) {
1178  RT fac = dyi * (face.isLow() ? RT(-1.0) : RT(1.));
1179  AMREX_HOST_DEVICE_PARALLEL_FOR_3D(b2d, i, j, k,
1180  {
1181  gp(i,j,k) = fac * (p(i,j,k) - p(i,j-1,k));
1182  });
1183  }
1184 #if (AMREX_SPACEDIM > 2)
1185  else {
1186  RT fac = dzi * (face.isLow() ? RT(-1.0) : RT(1.));
1187  AMREX_HOST_DEVICE_PARALLEL_FOR_3D(b2d, i, j, k,
1188  {
1189  gp(i,j,k) = fac * (p(i,j,k) - p(i,j,k-1));
1190  });
1191  }
1192 #endif
1193 #endif
1194  }
1195  }
1196  }
1197 }
1198 
1199 extern template class MLPoissonT<MultiFab>;
1200 
1202 
1203 }
1204 
1205 #endif
#define BL_PROFILE(a)
Definition: AMReX_BLProfiler.H:551
#define AMREX_ASSERT(EX)
Definition: AMReX_BLassert.H:38
#define AMREX_ALWAYS_ASSERT(EX)
Definition: AMReX_BLassert.H:50
#define AMREX_LAUNCH_HOST_DEVICE_LAMBDA(...)
Definition: AMReX_GpuLaunch.nolint.H:16
#define AMREX_HOST_DEVICE_PARALLEL_FOR_3D(...)
Definition: AMReX_GpuLaunch.nolint.H:54
#define AMREX_GPU_DEVICE
Definition: AMReX_GpuQualifiers.H:18
#define AMREX_D_TERM(a, b, c)
Definition: AMReX_SPACE.H:129
#define AMREX_D_DECL(a, b, c)
Definition: AMReX_SPACE.H:104
A collection of Boxes stored in an Array.
Definition: AMReX_BoxArray.H:549
Long size() const noexcept
Return the number of boxes in the BoxArray.
Definition: AMReX_BoxArray.H:596
AMREX_GPU_HOST_DEVICE IntVectND< dim > length() const noexcept
Return the length of the BoxND.
Definition: AMReX_Box.H:146
const Real * InvCellSize() const noexcept
Returns the inverse cellsize for each coordinate direction.
Definition: AMReX_CoordSys.H:82
Calculates the distribution of FABs to MPI processes.
Definition: AMReX_DistributionMapping.H:41
void define(const BoxArray &boxes, int nprocs=ParallelDescriptor::NProcs())
Build mapping out of BoxArray over nprocs processors. You need to call this if you built your Distrib...
Definition: AMReX_DistributionMapping.cpp:340
static DistributionMapping makeSFC(const MultiFab &weight, bool sort=true)
Definition: AMReX_DistributionMapping.cpp:1757
Definition: AMReX_FabFactory.H:50
Rectangular problem domain geometry.
Definition: AMReX_Geometry.H:73
Definition: AMReX_MFIter.H:57
Box tilebox() const noexcept
Return the tile Box at the current index.
Definition: AMReX_MFIter.cpp:385
bool isValid() const noexcept
Is the iterator valid i.e. is it associated with a FAB?
Definition: AMReX_MFIter.H:141
Definition: AMReX_MLALaplacian.H:14
Definition: AMReX_MLCellABecLap.H:13
typename MF::value_type RT
Definition: AMReX_MLCellABecLap.H:17
typename MLLinOpT< MF >::Location Location
Definition: AMReX_MLCellABecLap.H:19
typename MF::fab_type FAB
Definition: AMReX_MLCellABecLap.H:16
Vector< Vector< std::unique_ptr< iMultiFab > > > m_overset_mask
Definition: AMReX_MLCellABecLap.H:85
void define(const Vector< Geometry > &a_geom, const Vector< BoxArray > &a_grids, const Vector< DistributionMapping > &a_dmap, const LPInfo &a_info=LPInfo(), const Vector< FabFactory< FAB > const * > &a_factory={})
Definition: AMReX_MLCellABecLap.H:94
void prepareForSolve() override
Definition: AMReX_MLCellABecLap.H:247
BoxArray makeNGrids(int grid_size) const
Definition: AMReX_MLCellLinOp.H:897
LinOpBCType BCType
Definition: AMReX_MLCellLinOp.H:27
Vector< Vector< Array< MultiMask, 2 *AMREX_SPACEDIM > > > m_maskvals
Definition: AMReX_MLCellLinOp.H:200
Vector< Vector< BndryRegisterT< MF > > > m_undrrelxr
Definition: AMReX_MLCellLinOp.H:197
bool m_has_metric_term
Definition: AMReX_MLCellLinOp.H:146
bool m_use_gauss_seidel
Definition: AMReX_MLCellLinOp.H:206
T get_d0(T const &dx, T const &dy, T const &) const noexcept
Definition: AMReX_MLLinOp.H:714
bool doAgglomeration() const noexcept
Definition: AMReX_MLLinOp.H:669
Vector< Vector< BoxArray > > m_grids
Definition: AMReX_MLLinOp.H:599
Vector< Vector< DistributionMapping > > m_dmap
Definition: AMReX_MLLinOp.H:600
int verbose
Definition: AMReX_MLLinOp.H:577
IntVect m_coarse_data_crse_ratio
Definition: AMReX_MLLinOp.H:626
bool needsCoarseDataForBC() const noexcept
Needs coarse data for bc?
Definition: AMReX_MLLinOp.H:177
bool hasHiddenDimension() const noexcept
Definition: AMReX_MLLinOp.H:695
Vector< Array< BCType, AMREX_SPACEDIM > > m_lobc
Definition: AMReX_MLLinOp.H:562
int hiddenDirection() const noexcept
Definition: AMReX_MLLinOp.H:696
Vector< int > m_domain_covered
Definition: AMReX_MLLinOp.H:602
const MLLinOpT< MF > * m_parent
Definition: AMReX_MLLinOp.H:587
Vector< Vector< Geometry > > m_geom
first Vector is for amr level and second is mg level
Definition: AMReX_MLLinOp.H:598
Box compactify(Box const &b) const noexcept
Definition: AMReX_MLLinOp.H:1272
bool m_needs_coarse_data_for_bc
Definition: AMReX_MLLinOp.H:624
int maxorder
Definition: AMReX_MLLinOp.H:579
LPInfo info
Definition: AMReX_MLLinOp.H:575
Vector< Array< BCType, AMREX_SPACEDIM > > m_hibc
Definition: AMReX_MLLinOp.H:563
T get_d1(T const &, T const &dy, T const &dz) const noexcept
Definition: AMReX_MLLinOp.H:724
LinOpBCType m_coarse_fine_bc_type
Definition: AMReX_MLLinOp.H:625
int m_num_amr_levels
Definition: AMReX_MLLinOp.H:583
Definition: AMReX_MLPoisson.H:16
typename MF::value_type RT
Definition: AMReX_MLPoisson.H:20
MLPoissonT< MF > & operator=(const MLPoissonT< MF > &)=delete
MF const * getACoeffs(int, int) const final
Definition: AMReX_MLPoisson.H:70
void copyNSolveSolution(MF &dst, MF const &src) const final
Definition: AMReX_MLPoisson.H:1134
bool isBottomSingular() const final
Is the bottom of MG singular?
Definition: AMReX_MLPoisson.H:59
void prepareForSolve() final
Definition: AMReX_MLPoisson.H:140
void Fapply(int amrlev, int mglev, MF &out, const MF &in) const final
Definition: AMReX_MLPoisson.H:182
void normalize(int amrlev, int mglev, MF &mf) const final
Divide mf by the diagonal component of the operator. Used by bicgstab.
Definition: AMReX_MLPoisson.H:313
typename MF::fab_type FAB
Definition: AMReX_MLPoisson.H:19
bool isSingular(int amrlev) const final
Is it singular on given AMR level?
Definition: AMReX_MLPoisson.H:58
~MLPoissonT() override
MLPoissonT(const MLPoissonT< MF > &)=delete
bool supportNSolve() const final
Definition: AMReX_MLPoisson.H:1054
MLPoissonT(MLPoissonT< MF > &&)=delete
void define(const Vector< Geometry > &a_geom, const Vector< BoxArray > &a_grids, const Vector< DistributionMapping > &a_dmap, const LPInfo &a_info=LPInfo(), const Vector< FabFactory< FAB > const * > &a_factory={})
Definition: AMReX_MLPoisson.H:112
RT getAScalar() const final
Definition: AMReX_MLPoisson.H:68
RT getBScalar() const final
Definition: AMReX_MLPoisson.H:69
void get_dpdn_on_domain_faces(Array< MF *, AMREX_SPACEDIM > const &dpdn, MF const &phi)
Compute dphi/dn on domain faces after the solver has converged.
Definition: AMReX_MLPoisson.H:1141
MLPoissonT()=default
Vector< int > m_is_singular
Definition: AMReX_MLPoisson.H:86
void FFlux(int amrlev, const MFIter &mfi, const Array< FAB *, AMREX_SPACEDIM > &flux, const FAB &sol, Location loc, int face_only=0) const final
Definition: AMReX_MLPoisson.H:857
typename MLLinOpT< MF >::Location Location
Definition: AMReX_MLPoisson.H:23
Array< MF const *, AMREX_SPACEDIM > getBCoeffs(int, int) const final
Definition: AMReX_MLPoisson.H:71
void Fsmooth(int amrlev, int mglev, MF &sol, const MF &rhs, int redblack) const final
Definition: AMReX_MLPoisson.H:366
std::unique_ptr< MLLinOpT< MF > > makeNLinOp(int grid_size) const final
Definition: AMReX_MLPoisson.H:1065
Definition: AMReX_MultiMask.H:18
Array4< int const > array(const MFIter &mfi) const noexcept
Definition: AMReX_MultiMask.H:40
MultiArray4< int const > const_arrays() const noexcept
Definition: AMReX_MultiMask.H:48
An Iterator over the Orientation of Faces of a Box.
Definition: AMReX_Orientation.H:135
AMREX_GPU_HOST_DEVICE bool isValid() const noexcept
Is the iterator valid?
Definition: AMReX_Orientation.H:156
Encapsulation of the Orientation of the Faces of a Box.
Definition: AMReX_Orientation.H:29
AMREX_GPU_HOST_DEVICE int coordDir() const noexcept
Returns the coordinate direction.
Definition: AMReX_Orientation.H:83
AMREX_GPU_HOST_DEVICE bool isLow() const noexcept
Returns true if Orientation is low.
Definition: AMReX_Orientation.H:89
static const Periodicity & NonPeriodic() noexcept
Definition: AMReX_Periodicity.cpp:52
A Real vector in SpaceDim-dimensional space.
Definition: AMReX_RealVect.H:32
This class is a thin wrapper around std::vector. Unlike vector, Vector::operator[] provides bound che...
Definition: AMReX_Vector.H:27
void streamSynchronize() noexcept
Definition: AMReX_GpuDevice.H:237
bool inLaunchRegion() noexcept
Definition: AMReX_GpuControl.H:86
bool notInLaunchRegion() noexcept
Definition: AMReX_GpuControl.H:87
MPI_Comm CommunicatorSub() noexcept
sub-communicator for current frame
Definition: AMReX_ParallelContext.H:70
MPI_Comm Communicator() noexcept
Definition: AMReX_ParallelDescriptor.H:210
int NProcs() noexcept
return the number of MPI ranks local to the current Parallel Context
Definition: AMReX_ParallelDescriptor.H:243
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlpoisson_jacobi(int i, int j, int, Array4< T > const &phi, Array4< T const > const &rhs, Array4< T const > const &Ax, T dhx, T dhy, Array4< T const > const &f0, Array4< int const > const &m0, Array4< T const > const &f1, Array4< int const > const &m1, Array4< T const > const &f2, Array4< int const > const &m2, Array4< T const > const &f3, Array4< int const > const &m3, Box const &vbox) noexcept
Definition: AMReX_MLPoisson_2D_K.H:310
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlpoisson_adotx(int i, int j, Array4< T > const &y, Array4< T const > const &x, T dhx, T dhy) noexcept
Definition: AMReX_MLPoisson_2D_K.H:13
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlpoisson_flux_yface_m(Box const &box, Array4< T > const &fy, Array4< T const > const &sol, T dyinv, int ylen, T dx, T probxlo) noexcept
Definition: AMReX_MLPoisson_2D_K.H:174
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlpoisson_flux_y_m(Box const &box, Array4< T > const &fy, Array4< T const > const &sol, T dyinv, T dx, T probxlo) noexcept
Definition: AMReX_MLPoisson_2D_K.H:136
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlpoisson_gsrb(int i, int j, int, Array4< T > const &phi, Array4< T const > const &rhs, T dhx, T dhy, Array4< T const > const &f0, Array4< int const > const &m0, Array4< T const > const &f1, Array4< int const > const &m1, Array4< T const > const &f2, Array4< int const > const &m2, Array4< T const > const &f3, Array4< int const > const &m3, Box const &vbox, int redblack) noexcept
Definition: AMReX_MLPoisson_2D_K.H:197
@ max
Definition: AMReX_ParallelReduce.H:17
Definition: AMReX_Amr.cpp:49
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
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlpoisson_flux_yface(Box const &box, Array4< T > const &fy, Array4< T const > const &sol, T dyinv, int ylen) noexcept
Definition: AMReX_MLPoisson_3D_K.H:90
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlpoisson_jacobi(int i, int, int, Array4< T > const &phi, Array4< T const > const &rhs, Array4< T const > const &Ax, T dhx, Array4< T const > const &f0, Array4< int const > const &m0, Array4< T const > const &f1, Array4< int const > const &m1, Box const &vbox) noexcept
Definition: AMReX_MLPoisson_1D_K.H:193
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlpoisson_adotx_m(int i, Array4< T > const &y, Array4< T const > const &x, T dhx, T dx, T probxlo) noexcept
Definition: AMReX_MLPoisson_1D_K.H:32
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlpoisson_gsrb_m(int i, int, int, Array4< T > const &phi, Array4< T const > const &rhs, T dhx, Array4< T const > const &f0, Array4< int const > const &m0, Array4< T const > const &f1, Array4< int const > const &m1, Box const &vbox, int redblack, T dx, T probxlo) noexcept
Definition: AMReX_MLPoisson_1D_K.H:162
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlpoisson_flux_x(Box const &box, Array4< T > const &fx, Array4< T const > const &sol, T dxinv) noexcept
Definition: AMReX_MLPoisson_1D_K.H:43
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlpoisson_adotx(int i, Array4< T > const &y, Array4< T const > const &x, T dhx) noexcept
Definition: AMReX_MLPoisson_1D_K.H:9
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlpoisson_flux_y(Box const &box, Array4< T > const &fy, Array4< T const > const &sol, T dyinv) noexcept
Definition: AMReX_MLPoisson_3D_K.H:72
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlpoisson_flux_xface_m(Box const &box, Array4< T > const &fx, Array4< T const > const &sol, T dxinv, int xlen, T dx, T probxlo) noexcept
Definition: AMReX_MLPoisson_1D_K.H:86
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE BoxND< dim > bdryNode(const BoxND< dim > &b, Orientation face, int len=1) noexcept
Similar to bdryLo and bdryHi except that it operates on the given face of box b.
Definition: AMReX_Box.H:1549
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlpoisson_jacobi_os(int i, int, int, Array4< T > const &phi, Array4< T const > const &rhs, Array4< T const > const &Ax, Array4< int const > const &osm, T dhx, Array4< T const > const &f0, Array4< int const > const &m0, Array4< T const > const &f1, Array4< int const > const &m1, Box const &vbox) noexcept
Definition: AMReX_MLPoisson_1D_K.H:216
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlpoisson_flux_x_m(Box const &box, Array4< T > const &fx, Array4< T const > const &sol, T dxinv, T dx, T probxlo) noexcept
Definition: AMReX_MLPoisson_1D_K.H:57
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE BoxND< dim > bdryLo(const BoxND< dim > &b, int dir, int len=1) noexcept
Returns the edge-centered BoxND (in direction dir) defining the low side of BoxND b.
Definition: AMReX_Box.H:1502
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlpoisson_gsrb(int i, int, int, Array4< T > const &phi, Array4< T const > const &rhs, T dhx, Array4< T const > const &f0, Array4< int const > const &m0, Array4< T const > const &f1, Array4< int const > const &m1, Box const &vbox, int redblack) noexcept
Definition: AMReX_MLPoisson_1D_K.H:102
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Dim3 end(BoxND< dim > const &box) noexcept
Definition: AMReX_Box.H:1890
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlpoisson_flux_xface(Box const &box, Array4< T > const &fx, Array4< T const > const &sol, T dxinv, int xlen) noexcept
Definition: AMReX_MLPoisson_1D_K.H:73
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
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlpoisson_normalize(int i, int, int, Array4< T > const &x, T dhx, T dx, T probxlo) noexcept
Definition: AMReX_MLPoisson_1D_K.H:268
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE BoxND< dim > surroundingNodes(const BoxND< dim > &b, int dir) noexcept
Returns a BoxND with NODE based coordinates in direction dir that encloses BoxND b....
Definition: AMReX_Box.H:1399
bool TilingIfNotGPU() noexcept
Definition: AMReX_MFIter.H:12
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlpoisson_gsrb_os(int i, int, int, Array4< T > const &phi, Array4< T const > const &rhs, Array4< int const > const &osm, T dhx, Array4< T const > const &f0, Array4< int const > const &m0, Array4< T const > const &f1, Array4< int const > const &m1, Box const &vbox, int redblack) noexcept
Definition: AMReX_MLPoisson_1D_K.H:130
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Dim3 begin(BoxND< dim > const &box) noexcept
Definition: AMReX_Box.H:1881
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlpoisson_flux_zface(Box const &box, Array4< T > const &fz, Array4< T const > const &sol, T dzinv, int zlen) noexcept
Definition: AMReX_MLPoisson_3D_K.H:130
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlpoisson_adotx_os(int i, Array4< T > const &y, Array4< T const > const &x, Array4< int const > const &osm, T dhx) noexcept
Definition: AMReX_MLPoisson_1D_K.H:18
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlpoisson_flux_z(Box const &box, Array4< T > const &fz, Array4< T const > const &sol, T dzinv) noexcept
Definition: AMReX_MLPoisson_3D_K.H:112
std::array< T, N > Array
Definition: AMReX_Array.H:24
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlpoisson_jacobi_m(int i, int, int, Array4< T > const &phi, Array4< T const > const &rhs, Array4< T const > const &Ax, T dhx, Array4< T const > const &f0, Array4< int const > const &m0, Array4< T const > const &f1, Array4< int const > const &m1, Box const &vbox, T dx, T probxlo) noexcept
Definition: AMReX_MLPoisson_1D_K.H:242
parallel copy or add
Definition: AMReX_FabArrayBase.H:536
Definition: AMReX_MLLinOp.H:35
bool has_metric_term
Definition: AMReX_MLLinOp.H:43
Location
Definition: AMReX_MLLinOp.H:87
FabArray memory allocation information.
Definition: AMReX_FabArray.H:66
Definition: AMReX_MFIter.H:20
MFItInfo & EnableTiling(const IntVect &ts=FabArrayBase::mfiter_tile_size) noexcept
Definition: AMReX_MFIter.H:29
MFItInfo & SetDynamic(bool f) noexcept
Definition: AMReX_MFIter.H:34