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