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