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 AMREX_ALWAYS_ASSERT_WITH_MESSAGE(!this->m_has_metric_term,
204 "MLPoisson: overset mask does not support metric terms. Use MLABecLaplacian or LPInfo::setMetricTerm(false).");
205}
206
207template <typename MF>
208MLPoissonT<MF>::~MLPoissonT () = default;
209
210template <typename MF>
211void
212MLPoissonT<MF>::prepareForSolve ()
213{
214 BL_PROFILE("MLPoisson::prepareForSolve()");
215
216 MLCellABecLapT<MF>::prepareForSolve();
217
218 m_is_singular.clear();
219 m_is_singular.resize(this->m_num_amr_levels, false);
220 auto itlo = std::ranges::find(this->m_lobc[0], BCType::Dirichlet);
221 auto ithi = std::ranges::find(this->m_hibc[0], BCType::Dirichlet);
222 if (itlo == this->m_lobc[0].end() && ithi == this->m_hibc[0].end())
223 { // No Dirichlet
224 for (int alev = 0; alev < this->m_num_amr_levels; ++alev)
225 {
226 // For now this assumes that overset regions are treated as Dirichlet bc's
227 if (this->m_domain_covered[alev] && !this->m_overset_mask[alev][0])
228 {
229 m_is_singular[alev] = true;
230 }
231 }
232 }
233 if (!m_is_singular[0] && this->m_needs_coarse_data_for_bc &&
234 this->m_coarse_fine_bc_type == LinOpBCType::Neumann)
235 {
236 AMREX_ASSERT(this->m_overset_mask[0][0] == nullptr);
237 auto bbox = this->m_grids[0][0].minimalBox();
238 for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
239 if (this->m_lobc[0][idim] == LinOpBCType::Dirichlet) {
240 bbox.growLo(idim,1);
241 }
242 if (this->m_hibc[0][idim] == LinOpBCType::Dirichlet) {
243 bbox.growHi(idim,1);
244 }
245 }
246 if (this->m_geom[0][0].Domain().contains(bbox)) {
247 m_is_singular[0] = true;
248 }
249 }
250}
251
252template <typename MF>
253void
254MLPoissonT<MF>::Fapply (int amrlev, int mglev, MF& out, const MF& in) const
255{
256 BL_PROFILE("MLPoisson::Fapply()");
257
258 const Real* dxinv = this->m_geom[amrlev][mglev].InvCellSize();
259
260 AMREX_D_TERM(const RT dhx = RT(dxinv[0]*dxinv[0]);,
261 const RT dhy = RT(dxinv[1]*dxinv[1]);,
262 const RT dhz = RT(dxinv[2]*dxinv[2]););
263
264#if (AMREX_SPACEDIM == 3)
265 RT dh0 = this->get_d0(dhx, dhy, dhz);
266 RT dh1 = this->get_d1(dhx, dhy, dhz);
267#endif
268
269#if (AMREX_SPACEDIM < 3)
270 const RT dx = RT(this->m_geom[amrlev][mglev].CellSize(0));
271 const RT probxlo = RT(this->m_geom[amrlev][mglev].ProbLo(0));
272#endif
273
274#ifdef AMREX_USE_GPU
275 if (Gpu::inLaunchRegion() && out.isFusingCandidate() && !this->hasHiddenDimension()) {
276 auto const& xma = in.const_arrays();
277 auto const& yma = out.arrays();
278 if (this->m_overset_mask[amrlev][mglev]) {
279 AMREX_ASSERT(!this->m_has_metric_term);
280 const auto& osmma = this->m_overset_mask[amrlev][mglev]->const_arrays();
281 ParallelFor(out,
282 [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept
283 {
284 amrex::ignore_unused(j,k);
285 mlpoisson_adotx_os(AMREX_D_DECL(i,j,k), yma[box_no], xma[box_no], osmma[box_no],
286 AMREX_D_DECL(dhx,dhy,dhz));
287 });
288 } else {
289#if (AMREX_SPACEDIM < 3)
290 if (this->m_has_metric_term) {
291 ParallelFor(out,
292 [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept
293 {
294 amrex::ignore_unused(j,k);
295 mlpoisson_adotx_m(AMREX_D_DECL(i,j,k), yma[box_no], xma[box_no],
296 AMREX_D_DECL(dhx,dhy,dhz), dx, probxlo);
297 });
298 } else
299#endif
300 {
301 ParallelFor(out,
302 [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept
303 {
304 amrex::ignore_unused(j,k);
305 mlpoisson_adotx(AMREX_D_DECL(i,j,k), yma[box_no], xma[box_no],
306 AMREX_D_DECL(dhx,dhy,dhz));
307 });
308 }
309 }
310 if (!Gpu::inNoSyncRegion()) {
311 Gpu::streamSynchronize();
312 }
313 } else
314#endif
315 {
316#ifdef AMREX_USE_OMP
317#pragma omp parallel if (Gpu::notInLaunchRegion())
318#endif
319 for (MFIter mfi(out, TilingIfNotGPU()); mfi.isValid(); ++mfi)
320 {
321 const Box& bx = mfi.tilebox();
322 const auto& xfab = in.array(mfi);
323 const auto& yfab = out.array(mfi);
324
325 if (this->m_overset_mask[amrlev][mglev]) {
326 AMREX_ASSERT(!this->m_has_metric_term);
327 const auto& osm = this->m_overset_mask[amrlev][mglev]->const_array(mfi);
328 AMREX_HOST_DEVICE_PARALLEL_FOR_3D( bx, i, j, k,
329 {
330 amrex::ignore_unused(j,k);
331 mlpoisson_adotx_os(AMREX_D_DECL(i,j,k), yfab, xfab, osm,
332 AMREX_D_DECL(dhx,dhy,dhz));
333 });
334 } else {
335#if (AMREX_SPACEDIM == 3)
336 if (this->hasHiddenDimension()) {
337 Box const& bx2d = this->compactify(bx);
338 const auto& xfab2d = this->compactify(xfab);
339 const auto& yfab2d = this->compactify(yfab);
340 AMREX_HOST_DEVICE_PARALLEL_FOR_3D(bx2d, i, j, k,
341 {
342 amrex::ignore_unused(k);
343 TwoD::mlpoisson_adotx(i, j, yfab2d, xfab2d, dh0, dh1);
344 });
345 } else {
346 AMREX_HOST_DEVICE_PARALLEL_FOR_3D(bx, i, j, k,
347 {
348 mlpoisson_adotx(i, j, k, yfab, xfab, dhx, dhy, dhz);
349 });
350 }
351#elif (AMREX_SPACEDIM == 2)
352 if (this->m_has_metric_term) {
353 AMREX_HOST_DEVICE_PARALLEL_FOR_3D(bx, i, j, k,
354 {
355 amrex::ignore_unused(k);
356 mlpoisson_adotx_m(i, j, yfab, xfab, dhx, dhy, dx, probxlo);
357 });
358 } else {
359 AMREX_HOST_DEVICE_PARALLEL_FOR_3D(bx, i, j, k,
360 {
361 amrex::ignore_unused(k);
362 mlpoisson_adotx(i, j, yfab, xfab, dhx, dhy);
363 });
364 }
365#elif (AMREX_SPACEDIM == 1)
366 if (this->m_has_metric_term) {
367 AMREX_HOST_DEVICE_PARALLEL_FOR_3D(bx, i, j, k,
368 {
369 amrex::ignore_unused(j,k);
370 mlpoisson_adotx_m(i, yfab, xfab, dhx, dx, probxlo);
371 });
372 } else {
373 AMREX_HOST_DEVICE_PARALLEL_FOR_3D(bx, i, j, k,
374 {
375 amrex::ignore_unused(j,k);
376 mlpoisson_adotx(i, yfab, xfab, dhx);
377 });
378 }
379#endif
380 }
381 }
382 }
383}
384
385template <typename MF>
386void
387MLPoissonT<MF>::normalize (int amrlev, int mglev, MF& mf) const
388{
389 amrex::ignore_unused(amrlev,mglev,mf);
390#if (AMREX_SPACEDIM != 3)
391 BL_PROFILE("MLPoisson::normalize()");
392
393 if (!this->m_has_metric_term) { return; }
394
395 const Real* dxinv = this->m_geom[amrlev][mglev].InvCellSize();
396 AMREX_D_TERM(const RT dhx = RT(dxinv[0]*dxinv[0]);,
397 const RT dhy = RT(dxinv[1]*dxinv[1]);,
398 const RT dhz = RT(dxinv[2]*dxinv[2]););
399 const RT dx = RT(this->m_geom[amrlev][mglev].CellSize(0));
400 const RT probxlo = RT(this->m_geom[amrlev][mglev].ProbLo(0));
401
402#ifdef AMREX_USE_GPU
403 if (Gpu::inLaunchRegion() && mf.isFusingCandidate()) {
404 auto const& ma = mf.arrays();
405 ParallelFor(mf,
406 [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept
407 {
408 mlpoisson_normalize(i,j,k, ma[box_no], AMREX_D_DECL(dhx,dhy,dhz), dx, probxlo);
409 });
410 if (!Gpu::inNoSyncRegion()) {
411 Gpu::streamSynchronize();
412 }
413 } else
414#endif
415 {
416#ifdef AMREX_USE_OMP
417#pragma omp parallel if (Gpu::notInLaunchRegion())
418#endif
419 for (MFIter mfi(mf, TilingIfNotGPU()); mfi.isValid(); ++mfi)
420 {
421 const Box& bx = mfi.tilebox();
422 const auto& fab = mf.array(mfi);
423
424#if (AMREX_SPACEDIM == 2)
425 AMREX_HOST_DEVICE_PARALLEL_FOR_3D(bx, i, j, k,
426 {
427 mlpoisson_normalize(i,j,k, fab, dhx, dhy, dx, probxlo);
428 });
429#else
430 AMREX_HOST_DEVICE_PARALLEL_FOR_3D(bx, i, j, k,
431 {
432 mlpoisson_normalize(i,j,k, fab, dhx, dx, probxlo);
433 });
434#endif
435 }
436 }
437#endif
438}
439
440template <typename MF>
441void
442MLPoissonT<MF>::Fsmooth (int amrlev, int mglev, MF& sol, const MF& rhs, int redblack) const
443{
444 BL_PROFILE("MLPoisson::Fsmooth()");
445
446 MF Ax;
447 if (! this->m_use_gauss_seidel) { // jacobi
448 Ax.define(sol.boxArray(), sol.DistributionMap(), sol.nComp(), 0,
449 MFInfo().SetArena(The_Async_Arena()));
450 Fapply(amrlev, mglev, Ax, sol);
451 }
452
453 const auto& undrrelxr = this->m_undrrelxr[amrlev][mglev];
454 const auto& maskvals = this->m_maskvals [amrlev][mglev];
455
456 OrientationIter oitr;
457
458 const auto& f0 = undrrelxr[oitr()]; ++oitr;
459 const auto& f1 = undrrelxr[oitr()]; ++oitr;
460#if (AMREX_SPACEDIM > 1)
461 const auto& f2 = undrrelxr[oitr()]; ++oitr;
462 const auto& f3 = undrrelxr[oitr()]; ++oitr;
463#if (AMREX_SPACEDIM > 2)
464 const auto& f4 = undrrelxr[oitr()]; ++oitr;
465 const auto& f5 = undrrelxr[oitr()]; ++oitr;
466#endif
467#endif
468
469 const MultiMask& mm0 = maskvals[0];
470 const MultiMask& mm1 = maskvals[1];
471#if (AMREX_SPACEDIM > 1)
472 const MultiMask& mm2 = maskvals[2];
473 const MultiMask& mm3 = maskvals[3];
474#if (AMREX_SPACEDIM > 2)
475 const MultiMask& mm4 = maskvals[4];
476 const MultiMask& mm5 = maskvals[5];
477#endif
478#endif
479
480 const Real* dxinv = this->m_geom[amrlev][mglev].InvCellSize();
481 AMREX_D_TERM(const RT dhx = RT(dxinv[0]*dxinv[0]);,
482 const RT dhy = RT(dxinv[1]*dxinv[1]);,
483 const RT dhz = RT(dxinv[2]*dxinv[2]););
484
485#if (AMREX_SPACEDIM == 3)
486 RT dh0 = RT(this->get_d0(dhx, dhy, dhz));
487 RT dh1 = RT(this->get_d1(dhx, dhy, dhz));
488#endif
489
490#if (AMREX_SPACEDIM < 3)
491 const RT dx = RT(this->m_geom[amrlev][mglev].CellSize(0));
492 const RT probxlo = RT(this->m_geom[amrlev][mglev].ProbLo(0));
493#endif
494
495 MFItInfo mfi_info;
496 if (Gpu::notInLaunchRegion()) { mfi_info.EnableTiling().SetDynamic(true); }
497
498#ifdef AMREX_USE_GPU
499 // Gauss-Seidel uses a fused launch over the cells of one color.
500 if (Gpu::inLaunchRegion()
501 && (this->m_use_gauss_seidel || sol.isFusingCandidate())
502 && ! this->hasHiddenDimension())
503 {
504 const auto& m0ma = mm0.const_arrays();
505 const auto& m1ma = mm1.const_arrays();
506#if (AMREX_SPACEDIM > 1)
507 const auto& m2ma = mm2.const_arrays();
508 const auto& m3ma = mm3.const_arrays();
509#if (AMREX_SPACEDIM > 2)
510 const auto& m4ma = mm4.const_arrays();
511 const auto& m5ma = mm5.const_arrays();
512#endif
513#endif
514
515 const auto& solnma = sol.arrays();
516 const auto& rhsma = rhs.const_arrays();
517 const IntVect rhs_ng = rhs.nGrowVect();
518
519 const auto& f0ma = f0.const_arrays();
520 const auto& f1ma = f1.const_arrays();
521#if (AMREX_SPACEDIM > 1)
522 const auto& f2ma = f2.const_arrays();
523 const auto& f3ma = f3.const_arrays();
524#if (AMREX_SPACEDIM > 2)
525 const auto& f4ma = f4.const_arrays();
526 const auto& f5ma = f5.const_arrays();
527#endif
528#endif
529
530 if (this->m_overset_mask[amrlev][mglev]) {
531 AMREX_ASSERT(!this->m_has_metric_term);
532 const auto& osmma = this->m_overset_mask[amrlev][mglev]->const_arrays();
533 if (this->m_use_gauss_seidel) {
534 ParallelForRedBlack(sol, redblack,
535 [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept
536 {
537 Box vbx = amrex::grow(Box(rhsma[box_no]), -rhs_ng);
538 mlpoisson_gsrb_os(i, j, k, solnma[box_no], rhsma[box_no],
539 osmma[box_no], AMREX_D_DECL(dhx, dhy, dhz),
540 f0ma[box_no], m0ma[box_no],
541 f1ma[box_no], m1ma[box_no],
542#if (AMREX_SPACEDIM > 1)
543 f2ma[box_no], m2ma[box_no],
544 f3ma[box_no], m3ma[box_no],
545#if (AMREX_SPACEDIM > 2)
546 f4ma[box_no], m4ma[box_no],
547 f5ma[box_no], m5ma[box_no],
548#endif
549#endif
550 vbx, redblack);
551 });
552 } else {
553 const auto& axma = Ax.const_arrays();
554 ParallelFor(sol,
555 [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept
556 {
557 Box vbx = amrex::grow(Box(rhsma[box_no]), -rhs_ng);
558 mlpoisson_jacobi_os(i, j, k, solnma[box_no], rhsma[box_no],
559 axma[box_no], osmma[box_no],
560 AMREX_D_DECL(dhx, dhy, dhz),
561 f0ma[box_no], m0ma[box_no],
562 f1ma[box_no], m1ma[box_no],
563#if (AMREX_SPACEDIM > 1)
564 f2ma[box_no], m2ma[box_no],
565 f3ma[box_no], m3ma[box_no],
566#if (AMREX_SPACEDIM > 2)
567 f4ma[box_no], m4ma[box_no],
568 f5ma[box_no], m5ma[box_no],
569#endif
570#endif
571 vbx);
572 });
573 }
574 }
575#if (AMREX_SPACEDIM < 3)
576 else if (this->m_has_metric_term) {
577 if (this->m_use_gauss_seidel) {
578 ParallelForRedBlack(sol, redblack,
579 [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept
580 {
581 Box vbx = amrex::grow(Box(rhsma[box_no]), -rhs_ng);
582 mlpoisson_gsrb_m(i, j, k, solnma[box_no], rhsma[box_no],
583 AMREX_D_DECL(dhx, dhy, dhz),
584 f0ma[box_no], m0ma[box_no],
585 f1ma[box_no], m1ma[box_no],
586#if (AMREX_SPACEDIM > 1)
587 f2ma[box_no], m2ma[box_no],
588 f3ma[box_no], m3ma[box_no],
589#endif
590 vbx, redblack,
591 dx, probxlo);
592 });
593 } else {
594 const auto& axma = Ax.const_arrays();
595 ParallelFor(sol,
596 [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept
597 {
598 Box vbx = amrex::grow(Box(rhsma[box_no]), -rhs_ng);
599 mlpoisson_jacobi_m(i, j, k, solnma[box_no], rhsma[box_no],
600 axma[box_no], AMREX_D_DECL(dhx, dhy, dhz),
601 f0ma[box_no], m0ma[box_no],
602 f1ma[box_no], m1ma[box_no],
603#if (AMREX_SPACEDIM > 1)
604 f2ma[box_no], m2ma[box_no],
605 f3ma[box_no], m3ma[box_no],
606#endif
607 vbx, dx, probxlo);
608 });
609 }
610 }
611#endif
612 else {
613 if (this->m_use_gauss_seidel) {
614 ParallelForRedBlack(sol, redblack,
615 [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept
616 {
617 Box vbx = amrex::grow(Box(rhsma[box_no]), -rhs_ng);
618 mlpoisson_gsrb(i, j, k, solnma[box_no], rhsma[box_no],
619 AMREX_D_DECL(dhx, dhy, dhz),
620 f0ma[box_no], m0ma[box_no],
621 f1ma[box_no], m1ma[box_no],
622#if (AMREX_SPACEDIM > 1)
623 f2ma[box_no], m2ma[box_no],
624 f3ma[box_no], m3ma[box_no],
625#if (AMREX_SPACEDIM > 2)
626 f4ma[box_no], m4ma[box_no],
627 f5ma[box_no], m5ma[box_no],
628#endif
629#endif
630 vbx, redblack);
631 });
632 } else {
633 const auto& axma = Ax.const_arrays();
634 ParallelFor(sol,
635 [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept
636 {
637 Box vbx = amrex::grow(Box(rhsma[box_no]), -rhs_ng);
638 mlpoisson_jacobi(i, j, k, solnma[box_no], rhsma[box_no],
639 axma[box_no], AMREX_D_DECL(dhx, dhy, dhz),
640 f0ma[box_no], m0ma[box_no],
641 f1ma[box_no], m1ma[box_no],
642#if (AMREX_SPACEDIM > 1)
643 f2ma[box_no], m2ma[box_no],
644 f3ma[box_no], m3ma[box_no],
645#if (AMREX_SPACEDIM > 2)
646 f4ma[box_no], m4ma[box_no],
647 f5ma[box_no], m5ma[box_no],
648#endif
649#endif
650 vbx);
651 });
652 }
653 }
654 if (!Gpu::inNoSyncRegion()) {
655 Gpu::streamSynchronize();
656 }
657 } else
658#endif
659 {
660#ifdef AMREX_USE_OMP
661#pragma omp parallel if (Gpu::notInLaunchRegion())
662#endif
663 for (MFIter mfi(sol,mfi_info); mfi.isValid(); ++mfi)
664 {
665 const auto& m0 = mm0.array(mfi);
666 const auto& m1 = mm1.array(mfi);
667#if (AMREX_SPACEDIM > 1)
668 const auto& m2 = mm2.array(mfi);
669 const auto& m3 = mm3.array(mfi);
670#if (AMREX_SPACEDIM > 2)
671 const auto& m4 = mm4.array(mfi);
672 const auto& m5 = mm5.array(mfi);
673#endif
674#endif
675
676 const Box& tbx = mfi.tilebox();
677 const Box& vbx = mfi.validbox();
678 const auto& solnfab = sol.array(mfi);
679 const auto& rhsfab = rhs.array(mfi);
680
681 const auto& f0fab = f0.array(mfi);
682 const auto& f1fab = f1.array(mfi);
683#if (AMREX_SPACEDIM > 1)
684 const auto& f2fab = f2.array(mfi);
685 const auto& f3fab = f3.array(mfi);
686#if (AMREX_SPACEDIM > 2)
687 const auto& f4fab = f4.array(mfi);
688 const auto& f5fab = f5.array(mfi);
689#endif
690#endif
691
692#if (AMREX_SPACEDIM == 1)
693 if (this->m_overset_mask[amrlev][mglev]) {
694 AMREX_ASSERT(!this->m_has_metric_term);
695 const auto& osm = this->m_overset_mask[amrlev][mglev]->const_array(mfi);
696 if (this->m_use_gauss_seidel) {
697 AMREX_HOST_DEVICE_PARALLEL_FOR_3D ( tbx, i, j, k,
698 {
699 mlpoisson_gsrb_os(i, j, k, solnfab, rhsfab, osm, dhx,
700 f0fab, m0,
701 f1fab, m1,
702 vbx, redblack);
703 });
704 } else {
705 const auto& axfab = Ax.const_array(mfi);
706 AMREX_HOST_DEVICE_PARALLEL_FOR_3D ( tbx, i, j, k,
707 {
708 mlpoisson_jacobi_os(i, j, k, solnfab, rhsfab, axfab,
709 osm, dhx,
710 f0fab, m0,
711 f1fab, m1,
712 vbx);
713 });
714 }
715 } else if (this->m_has_metric_term) {
716 if (this->m_use_gauss_seidel) {
717 AMREX_HOST_DEVICE_PARALLEL_FOR_3D ( tbx, i, j, k,
718 {
719 mlpoisson_gsrb_m(i, j, k, solnfab, rhsfab, dhx,
720 f0fab, m0,
721 f1fab, m1,
722 vbx, redblack,
723 dx, probxlo);
724 });
725 } else {
726 const auto& axfab = Ax.const_array(mfi);
727 AMREX_HOST_DEVICE_PARALLEL_FOR_3D ( tbx, i, j, k,
728 {
729 mlpoisson_jacobi_m(i, j, k, solnfab, rhsfab, axfab, dhx,
730 f0fab, m0,
731 f1fab, m1,
732 vbx, dx, probxlo);
733 });
734 }
735 } else {
736 if (this->m_use_gauss_seidel) {
737 AMREX_HOST_DEVICE_PARALLEL_FOR_3D ( tbx, i, j, k,
738 {
739 mlpoisson_gsrb(i, j, k, solnfab, rhsfab, dhx,
740 f0fab, m0,
741 f1fab, m1,
742 vbx, redblack);
743 });
744 } else {
745 const auto& axfab = Ax.const_array(mfi);
746 AMREX_HOST_DEVICE_PARALLEL_FOR_3D ( tbx, i, j, k,
747 {
748 mlpoisson_jacobi(i, j, k, solnfab, rhsfab, axfab, dhx,
749 f0fab, m0,
750 f1fab, m1,
751 vbx);
752 });
753 }
754 }
755#endif
756
757#if (AMREX_SPACEDIM == 2)
758 if (this->m_overset_mask[amrlev][mglev]) {
759 AMREX_ASSERT(!this->m_has_metric_term);
760 const auto& osm = this->m_overset_mask[amrlev][mglev]->const_array(mfi);
761 if (this->m_use_gauss_seidel) {
762 AMREX_HOST_DEVICE_PARALLEL_FOR_3D ( tbx, i, j, k,
763 {
764 mlpoisson_gsrb_os(i, j, k, solnfab, rhsfab, osm, dhx, dhy,
765 f0fab, m0,
766 f1fab, m1,
767 f2fab, m2,
768 f3fab, m3,
769 vbx, redblack);
770 });
771 } else {
772 const auto& axfab = Ax.const_array(mfi);
773 AMREX_HOST_DEVICE_PARALLEL_FOR_3D ( tbx, i, j, k,
774 {
775 mlpoisson_jacobi_os(i, j, k, solnfab, rhsfab, axfab,
776 osm, dhx, dhy,
777 f0fab, m0,
778 f1fab, m1,
779 f2fab, m2,
780 f3fab, m3,
781 vbx);
782 });
783 }
784 } else if (this->m_has_metric_term) {
785 if (this->m_use_gauss_seidel) {
786 AMREX_HOST_DEVICE_PARALLEL_FOR_3D ( tbx, i, j, k,
787 {
788 mlpoisson_gsrb_m(i, j, k, solnfab, rhsfab, dhx, dhy,
789 f0fab, m0,
790 f1fab, m1,
791 f2fab, m2,
792 f3fab, m3,
793 vbx, redblack,
794 dx, probxlo);
795 });
796 } else {
797 const auto& axfab = Ax.const_array(mfi);
798 AMREX_HOST_DEVICE_PARALLEL_FOR_3D ( tbx, i, j, k,
799 {
800 mlpoisson_jacobi_m(i, j, k, solnfab, rhsfab, axfab, dhx, dhy,
801 f0fab, m0,
802 f1fab, m1,
803 f2fab, m2,
804 f3fab, m3,
805 vbx, dx, probxlo);
806 });
807 }
808 } else {
809 if (this->m_use_gauss_seidel) {
810 AMREX_HOST_DEVICE_PARALLEL_FOR_3D ( tbx, i, j, k,
811 {
812 mlpoisson_gsrb(i, j, k, solnfab, rhsfab, dhx, dhy,
813 f0fab, m0,
814 f1fab, m1,
815 f2fab, m2,
816 f3fab, m3,
817 vbx, redblack);
818 });
819 } else {
820 const auto& axfab = Ax.const_array(mfi);
821 AMREX_HOST_DEVICE_PARALLEL_FOR_3D ( tbx, i, j, k,
822 {
823 mlpoisson_jacobi(i, j, k, solnfab, rhsfab, axfab, dhx, dhy,
824 f0fab, m0,
825 f1fab, m1,
826 f2fab, m2,
827 f3fab, m3,
828 vbx);
829 });
830 }
831 }
832#endif
833
834#if (AMREX_SPACEDIM == 3)
835 if (this->m_overset_mask[amrlev][mglev]) {
836 AMREX_ASSERT(!this->m_has_metric_term);
837 const auto& osm = this->m_overset_mask[amrlev][mglev]->const_array(mfi);
838 if (this->m_use_gauss_seidel) {
839 AMREX_HOST_DEVICE_PARALLEL_FOR_3D ( tbx, i, j, k,
840 {
841 mlpoisson_gsrb_os(i, j, k, solnfab, rhsfab, osm, dhx, dhy, dhz,
842 f0fab, m0,
843 f1fab, m1,
844 f2fab, m2,
845 f3fab, m3,
846 f4fab, m4,
847 f5fab, m5,
848 vbx, redblack);
849 });
850 } else {
851 const auto& axfab = Ax.const_array(mfi);
852 AMREX_HOST_DEVICE_PARALLEL_FOR_3D ( tbx, i, j, k,
853 {
854 mlpoisson_jacobi_os(i, j, k, solnfab, rhsfab, axfab,
855 osm, dhx, dhy, dhz,
856 f0fab, m0,
857 f1fab, m1,
858 f2fab, m2,
859 f3fab, m3,
860 f4fab, m4,
861 f5fab, m5,
862 vbx);
863 });
864 }
865 } else if (this->hasHiddenDimension()) {
866 Box const& tbx_2d = this->compactify(tbx);
867 Box const& vbx_2d = this->compactify(vbx);
868 const auto& solnfab_2d = this->compactify(solnfab);
869 const auto& rhsfab_2d = this->compactify(rhsfab);
870 const auto& f0fab_2d = this->compactify(this->get_d0(f0fab,f1fab,f2fab));
871 const auto& f1fab_2d = this->compactify(this->get_d1(f0fab,f1fab,f2fab));
872 const auto& f2fab_2d = this->compactify(this->get_d0(f3fab,f4fab,f5fab));
873 const auto& f3fab_2d = this->compactify(this->get_d1(f3fab,f4fab,f5fab));
874 const auto& m0_2d = this->compactify(this->get_d0(m0,m1,m2));
875 const auto& m1_2d = this->compactify(this->get_d1(m0,m1,m2));
876 const auto& m2_2d = this->compactify(this->get_d0(m3,m4,m5));
877 const auto& m3_2d = this->compactify(this->get_d1(m3,m4,m5));
878 if (this->m_use_gauss_seidel) {
879 AMREX_HOST_DEVICE_PARALLEL_FOR_3D ( tbx_2d, i, j, k,
880 {
881 TwoD::mlpoisson_gsrb(i, j, k, solnfab_2d, rhsfab_2d, dh0, dh1,
882 f0fab_2d, m0_2d,
883 f1fab_2d, m1_2d,
884 f2fab_2d, m2_2d,
885 f3fab_2d, m3_2d,
886 vbx_2d, redblack);
887 });
888 } else {
889 const auto& axfab = Ax.const_array(mfi);
890 const auto& axfab_2d = this->compactify(axfab);
891 AMREX_HOST_DEVICE_PARALLEL_FOR_3D ( tbx_2d, i, j, k,
892 {
893 TwoD::mlpoisson_jacobi(i, j, k, solnfab_2d, rhsfab_2d,
894 axfab_2d, dh0, dh1,
895 f0fab_2d, m0_2d,
896 f1fab_2d, m1_2d,
897 f2fab_2d, m2_2d,
898 f3fab_2d, m3_2d,
899 vbx_2d);
900 });
901 }
902 } else {
903 if (this->m_use_gauss_seidel) {
904 AMREX_HOST_DEVICE_PARALLEL_FOR_3D ( tbx, i, j, k,
905 {
906 mlpoisson_gsrb(i, j, k, solnfab, rhsfab, dhx, dhy, dhz,
907 f0fab, m0,
908 f1fab, m1,
909 f2fab, m2,
910 f3fab, m3,
911 f4fab, m4,
912 f5fab, m5,
913 vbx, redblack);
914 });
915 } else {
916 const auto& axfab = Ax.const_array(mfi);
917 AMREX_HOST_DEVICE_PARALLEL_FOR_3D ( tbx, i, j, k,
918 {
919 mlpoisson_jacobi(i, j, k, solnfab, rhsfab, axfab,
920 dhx, dhy, dhz,
921 f0fab, m0,
922 f1fab, m1,
923 f2fab, m2,
924 f3fab, m3,
925 f4fab, m4,
926 f5fab, m5,
927 vbx);
928 });
929 }
930 }
931#endif
932 }
933 }
934}
935
936template <typename MF>
937void
938MLPoissonT<MF>::FFlux (int amrlev, const MFIter& mfi,
939 const Array<FAB*,AMREX_SPACEDIM>& flux,
940 const FAB& sol, Location, const int face_only) const
941{
942 BL_PROFILE("MLPoisson::FFlux()");
943
944 const int mglev = 0;
945 const Box& box = mfi.tilebox();
946 const Real* dxinv = this->m_geom[amrlev][mglev].InvCellSize();
947
948 AMREX_D_TERM(const auto& fxarr = flux[0]->array();,
949 const auto& fyarr = flux[1]->array();,
950 const auto& fzarr = flux[2]->array(););
951 const auto& solarr = sol.array();
952
953#if (AMREX_SPACEDIM != 3)
954 const RT dx = RT(this->m_geom[amrlev][mglev].CellSize(0));
955 const RT probxlo = RT(this->m_geom[amrlev][mglev].ProbLo(0));
956#endif
957
958#if (AMREX_SPACEDIM == 3)
959 if (face_only) {
960 if (this->hiddenDirection() != 0) {
961 RT fac = RT(dxinv[0]);
962 Box blo = amrex::bdryLo(box, 0);
963 int blen = box.length(0);
964 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( blo, tbox,
965 {
966 mlpoisson_flux_xface(tbox, fxarr, solarr, fac, blen);
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 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( blo, tbox,
976 {
977 mlpoisson_flux_yface(tbox, fyarr, solarr, fac, blen);
978 });
979 } else {
980 flux[1]->template setVal<RunOn::Device>(RT(0.0));
981 }
982 if (this->hiddenDirection() != 2) {
983 RT fac = RT(dxinv[2]);
984 Box blo = amrex::bdryLo(box, 2);
985 int blen = box.length(2);
986 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( blo, tbox,
987 {
988 mlpoisson_flux_zface(tbox, fzarr, solarr, fac, blen);
989 });
990 } else {
991 flux[2]->template setVal<RunOn::Device>(RT(0.0));
992 }
993 } else {
994 if (this->hiddenDirection() != 0) {
995 RT fac = RT(dxinv[0]);
996 Box bflux = amrex::surroundingNodes(box, 0);
997 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( bflux, tbox,
998 {
999 mlpoisson_flux_x(tbox, fxarr, solarr, fac);
1000 });
1001 } else {
1002 flux[0]->template setVal<RunOn::Device>(RT(0.0));
1003 }
1004 if (this->hiddenDirection() != 1) {
1005 RT fac = RT(dxinv[1]);
1006 Box bflux = amrex::surroundingNodes(box, 1);
1007 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( bflux, tbox,
1008 {
1009 mlpoisson_flux_y(tbox, fyarr, solarr, fac);
1010 });
1011 } else {
1012 flux[1]->template setVal<RunOn::Device>(RT(0.0));
1013 }
1014 if (this->hiddenDirection() != 2) {
1015 RT fac = RT(dxinv[2]);
1016 Box bflux = amrex::surroundingNodes(box, 2);
1017 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( bflux, tbox,
1018 {
1019 mlpoisson_flux_z(tbox, fzarr, solarr, fac);
1020 });
1021 } else {
1022 flux[2]->template setVal<RunOn::Device>(RT(0.0));
1023 }
1024 }
1025#elif (AMREX_SPACEDIM == 2)
1026 if (face_only) {
1027 if (this->hiddenDirection() != 0) {
1028 RT fac = RT(dxinv[0]);
1029 Box blo = amrex::bdryLo(box, 0);
1030 int blen = box.length(0);
1031 if (this->m_has_metric_term) {
1032 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( blo, tbox,
1033 {
1034 mlpoisson_flux_xface_m(tbox, fxarr, solarr, fac, blen, dx, probxlo);
1035 });
1036 } else {
1037 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( blo, tbox,
1038 {
1039 mlpoisson_flux_xface(tbox, fxarr, solarr, fac, blen);
1040 });
1041 }
1042 } else {
1043 flux[0]->template setVal<RunOn::Device>(RT(0.0));
1044 }
1045 if (this->hiddenDirection() != 1) {
1046 RT fac = RT(dxinv[1]);
1047 Box blo = amrex::bdryLo(box, 1);
1048 int blen = box.length(1);
1049 if (this->m_has_metric_term) {
1050 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( blo, tbox,
1051 {
1052 mlpoisson_flux_yface_m(tbox, fyarr, solarr, fac, blen, dx, probxlo);
1053 });
1054 } else {
1055 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( blo, tbox,
1056 {
1057 mlpoisson_flux_yface(tbox, fyarr, solarr, fac, blen);
1058 });
1059 }
1060 } else {
1061 flux[1]->template setVal<RunOn::Device>(RT(0.0));
1062 }
1063 } else {
1064 if (this->hiddenDirection() != 0) {
1065 RT fac = RT(dxinv[0]);
1066 Box bflux = amrex::surroundingNodes(box, 0);
1067 if (this->m_has_metric_term) {
1068 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( bflux, tbox,
1069 {
1070 mlpoisson_flux_x_m(tbox, fxarr, solarr, fac, dx, probxlo);
1071 });
1072 } else {
1073 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( bflux, tbox,
1074 {
1075 mlpoisson_flux_x(tbox, fxarr, solarr, fac);
1076 });
1077 }
1078 } else {
1079 flux[0]->template setVal<RunOn::Device>(RT(0.0));
1080 }
1081 if (this->hiddenDirection() != 1) {
1082 RT fac = RT(dxinv[1]);
1083 Box bflux = amrex::surroundingNodes(box, 1);
1084 if (this->m_has_metric_term) {
1085 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( bflux, tbox,
1086 {
1087 mlpoisson_flux_y_m(tbox, fyarr, solarr, fac, dx, probxlo);
1088 });
1089 } else {
1090 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( bflux, tbox,
1091 {
1092 mlpoisson_flux_y(tbox, fyarr, solarr, fac);
1093 });
1094 }
1095 } else {
1096 flux[1]->template setVal<RunOn::Device>(RT(0.0));
1097 }
1098 }
1099#else
1100 if (face_only) {
1101 RT fac = RT(dxinv[0]);
1102 Box blo = amrex::bdryLo(box, 0);
1103 int blen = box.length(0);
1104 if (this->m_has_metric_term) {
1105 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( blo, tbox,
1106 {
1107 mlpoisson_flux_xface_m(tbox, fxarr, solarr, fac, blen, dx, probxlo);
1108 });
1109 } else {
1110 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( blo, tbox,
1111 {
1112 mlpoisson_flux_xface(tbox, fxarr, solarr, fac, blen);
1113 });
1114 }
1115 } else {
1116 RT fac = RT(dxinv[0]);
1117 Box bflux = amrex::surroundingNodes(box, 0);
1118 if (this->m_has_metric_term) {
1119 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( bflux, tbox,
1120 {
1121 mlpoisson_flux_x_m(tbox, fxarr, solarr, fac, dx, probxlo);
1122 });
1123 } else {
1124 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( bflux, tbox,
1125 {
1126 mlpoisson_flux_x(tbox, fxarr, solarr, fac);
1127 });
1128 }
1129 }
1130#endif
1131}
1132
1133template <typename MF>
1134bool
1135MLPoissonT<MF>::supportNSolve () const
1136{
1137 bool support = true;
1138 if (this->m_domain_covered[0]) { support = false; }
1139 if (this->doAgglomeration()) { support = false; }
1140 if (AMREX_SPACEDIM != 3) { support = false; }
1141 return support;
1142}
1143
1144template <typename MF>
1145std::unique_ptr<MLLinOpT<MF>>
1146MLPoissonT<MF>::makeNLinOp (int grid_size) const
1147{
1148 const Geometry& geom = this->m_geom[0].back();
1149 const BoxArray& ba = this->makeNGrids(grid_size);
1150
1151 DistributionMapping dm;
1152 {
1153 const std::vector<std::vector<int> >& sfc = DistributionMapping::makeSFC(ba);
1154 Vector<int> pmap(ba.size());
1155 AMREX_ALWAYS_ASSERT(ParallelContext::CommunicatorSub() == ParallelDescriptor::Communicator());
1156 const int nprocs = ParallelDescriptor::NProcs();
1157 for (int iproc = 0; iproc < nprocs; ++iproc) {
1158 for (int ibox : sfc[iproc]) {
1159 pmap[ibox] = iproc;
1160 }
1161 }
1162 dm.define(std::move(pmap));
1163 }
1164
1165 LPInfo minfo{};
1166 minfo.has_metric_term = this->info.has_metric_term;
1167
1168 std::unique_ptr<MLLinOpT<MF>> r{new MLALaplacianT<MF>({geom}, {ba}, {dm}, minfo)};
1169 auto nop = dynamic_cast<MLALaplacianT<MF>*>(r.get());
1170 if (!nop) {
1171 return nullptr;
1172 }
1173
1174 nop->m_parent = this;
1175
1176 nop->setMaxOrder(this->maxorder);
1177 nop->setVerbose(this->verbose);
1178
1179 nop->setDomainBC(this->m_lobc, this->m_hibc);
1180
1181 if (this->needsCoarseDataForBC())
1182 {
1183 const Real* dx0 = this->m_geom[0][0].CellSize();
1184 RealVect fac(this->m_coarse_data_crse_ratio);
1185 fac *= Real(0.5);
1186 RealVect cbloc {AMREX_D_DECL(dx0[0]*fac[0], dx0[1]*fac[1], dx0[2]*fac[2])};
1187 nop->setCoarseFineBCLocation(cbloc);
1188 }
1189
1190 nop->setScalars(1.0, -1.0);
1191
1192 const Real* dxinv = geom.InvCellSize();
1193 RT dxscale = RT(dxinv[0]);
1194#if (AMREX_SPACEDIM >= 2)
1195 dxscale = std::max(dxscale,RT(dxinv[1]));
1196#endif
1197#if (AMREX_SPACEDIM == 3)
1198 dxscale = std::max(dxscale,RT(dxinv[2]));
1199#endif
1200
1201 MF alpha(ba, dm, 1, 0, MFInfo().SetArena(The_Async_Arena()));
1202 alpha.setVal(RT(1.e30)*dxscale*dxscale);
1203
1204 MF foo(this->m_grids[0].back(), this->m_dmap[0].back(), 1, 0, MFInfo().SetAlloc(false));
1205 const FabArrayBase::CPC& cpc = alpha.getCPC(IntVect(0),foo,IntVect(0),Periodicity::NonPeriodic());
1206 alpha.setVal(RT(0.0), cpc, 0, 1);
1207
1208 nop->setACoeffs(0, alpha);
1209
1210 return r;
1211}
1212
1213template <typename MF>
1214void
1215MLPoissonT<MF>::copyNSolveSolution (MF& dst, MF const& src) const
1216{
1217 dst.ParallelCopy(src);
1218}
1219
1220template <typename MF>
1221void
1222MLPoissonT<MF>::get_dpdn_on_domain_faces (Array<MF*,AMREX_SPACEDIM> const& dpdn,
1223 MF const& phi)
1224{
1225 BL_PROFILE("MLPoisson::dpdn_faces()");
1226
1227 // We do not need to call applyBC because this function is used by the
1228 // OpenBC solver after solver has converged. That means the BC has been
1229 // filled to check the residual.
1230
1231 Box const& domain0 = this->m_geom[0][0].Domain();
1232 AMREX_D_TERM(const RT dxi = RT(this->m_geom[0][0].InvCellSize(0));,
1233 const RT dyi = RT(this->m_geom[0][0].InvCellSize(1));,
1234 const RT dzi = RT(this->m_geom[0][0].InvCellSize(2));)
1235
1236#ifdef AMREX_USE_OMP
1237#pragma omp parallel if (Gpu::notInLaunchRegion())
1238#endif
1239 for (MFIter mfi(phi); mfi.isValid(); ++mfi)
1240 {
1241 Box const& vbx = mfi.validbox();
1242 for (OrientationIter oit; oit.isValid(); ++oit) {
1243 Orientation face = oit();
1244 if (vbx[face] == domain0[face]) {
1245 int dir = face.coordDir();
1246 auto const& p = phi.const_array(mfi);
1247 auto const& gp = dpdn[dir]->array(mfi);
1248 Box const& b2d = amrex::bdryNode(vbx,face);
1249 if (dir == 0) {
1250 // because it's dphi/dn, not dphi/dx.
1251 RT fac = dxi * (face.isLow() ? RT(-1.0) : RT(1.));
1252 AMREX_HOST_DEVICE_PARALLEL_FOR_3D(b2d, i, j, k,
1253 {
1254 gp(i,j,k) = fac * (p(i,j,k) - p(i-1,j,k));
1255 });
1256 }
1257#if (AMREX_SPACEDIM > 1)
1258 else if (dir == 1) {
1259 RT fac = dyi * (face.isLow() ? RT(-1.0) : RT(1.));
1260 AMREX_HOST_DEVICE_PARALLEL_FOR_3D(b2d, i, j, k,
1261 {
1262 gp(i,j,k) = fac * (p(i,j,k) - p(i,j-1,k));
1263 });
1264 }
1265#if (AMREX_SPACEDIM > 2)
1266 else {
1267 RT fac = dzi * (face.isLow() ? RT(-1.0) : RT(1.));
1268 AMREX_HOST_DEVICE_PARALLEL_FOR_3D(b2d, i, j, k,
1269 {
1270 gp(i,j,k) = fac * (p(i,j,k) - p(i,j,k-1));
1271 });
1272 }
1273#endif
1274#endif
1275 }
1276 }
1277 }
1278}
1279
1280extern template class MLPoissonT<MultiFab>;
1281
1283using MLPoisson = MLPoissonT<MultiFab>;
1284
1285}
1286
1287#endif
BL_PROFILE
#define BL_PROFILE(a)
Definition AMReX_BLProfiler.H:562
AMREX_ALWAYS_ASSERT_WITH_MESSAGE
#define AMREX_ALWAYS_ASSERT_WITH_MESSAGE(EX, MSG)
Definition AMReX_BLassert.H:49
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:681
amrex::BoxArray::size
Long size() const noexcept
Return the number of boxes in the BoxArray.
Definition AMReX_BoxArray.H:758
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:1770
amrex::FabFactory
Abstract factory interface for creating, aliasing, and destroying FAB objects.
Definition AMReX_FabFactory.H:73
amrex::Geometry
Rectangular problem domain geometry.
Definition AMReX_Geometry.H:85
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:176
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:151
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:160
amrex::MLCellABecLapT::prepareForSolve
void prepareForSolve() override
Standard hook called before MLMG iterates (fixes BC data, etc.).
Definition AMReX_MLCellABecLap.H:319
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:1206
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:1106
amrex::MLLinOpT::doAgglomeration
bool doAgglomeration() const noexcept
Definition AMReX_MLLinOp.H:1038
amrex::MLLinOpT::m_hibc
Vector< Array< BCType, 3 > > m_hibc
Definition AMReX_MLLinOp.H:915
amrex::MLLinOpT::m_grids
Vector< Vector< BoxArray > > m_grids
Definition AMReX_MLLinOp.H:957
amrex::MLLinOpT::m_dmap
Vector< Vector< DistributionMapping > > m_dmap
Definition AMReX_MLLinOp.H:958
amrex::MLLinOpT::verbose
int verbose
Definition AMReX_MLLinOp.H:932
amrex::MLLinOpT::m_coarse_data_crse_ratio
IntVect m_coarse_data_crse_ratio
Definition AMReX_MLLinOp.H:986
amrex::MLLinOpT::needsCoarseDataForBC
bool needsCoarseDataForBC() const noexcept
Needs coarse data for bc?
Definition AMReX_MLLinOp.H:230
amrex::MLLinOpT::hasHiddenDimension
bool hasHiddenDimension() const noexcept
Definition AMReX_MLLinOp.H:1087
amrex::MLLinOpT::hiddenDirection
int hiddenDirection() const noexcept
Definition AMReX_MLLinOp.H:1088
amrex::MLLinOpT::m_lobc
Vector< Array< BCType, 3 > > m_lobc
Definition AMReX_MLLinOp.H:914
amrex::MLLinOpT::m_domain_covered
Vector< int > m_domain_covered
Definition AMReX_MLLinOp.H:960
amrex::MLLinOpT::m_parent
const MLLinOpT< MF > * m_parent
Definition AMReX_MLLinOp.H:945
amrex::MLLinOpT::m_geom
Vector< Vector< Geometry > > m_geom
first Vector is for amr level and second is mg level
Definition AMReX_MLLinOp.H:956
amrex::MLLinOpT::compactify
Box compactify(Box const &b) const noexcept
Definition AMReX_MLLinOp.H:1842
amrex::MLLinOpT::m_needs_coarse_data_for_bc
bool m_needs_coarse_data_for_bc
Definition AMReX_MLLinOp.H:984
amrex::MLLinOpT::maxorder
int maxorder
Definition AMReX_MLLinOp.H:935
amrex::MLLinOpT::info
LPInfo info
Definition AMReX_MLLinOp.H:930
amrex::MLLinOpT::get_d1
T get_d1(T const &, T const &dy, T const &dz) const noexcept
Definition AMReX_MLLinOp.H:1116
amrex::MLLinOpT::m_coarse_fine_bc_type
LinOpBCType m_coarse_fine_bc_type
Definition AMReX_MLLinOp.H:985
amrex::MLLinOpT::m_num_amr_levels
int m_num_amr_levels
Definition AMReX_MLLinOp.H:941
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:1215
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:1222
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:212
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:254
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:387
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:1135
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:938
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:442
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:1146
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:80
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::grow
__host__ __device__ BoxND< dim > grow(const BoxND< dim > &b, int i) noexcept
Return a copy of b grown uniformly by i cells in every direction.
Definition AMReX_Box.H:1326
amrex::Array
std::array< T, N > Array
Definition AMReX_Array.H:31
amrex::The_Async_Arena
Arena * The_Async_Arena()
Definition AMReX_Arena.cpp:839
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:273
amrex::ParallelForRedBlack
void ParallelForRedBlack(MF const &mf, int redblack, F &&f)
ParallelFor over the red or black points of the valid region of a MultiFab/FabArray.
Definition AMReX_MFParallelFor.H:529
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::Box
BoxND< 3 > Box
Box is an alias for amrex::BoxND instantiated with AMREX_SPACEDIM.
Definition AMReX_BaseFwd.H:35
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:630
amrex::LPInfo
Configuration knobs for multilevel linear operators (grid agglomeration, metrics, etc....
Definition AMReX_MLLinOp.H:53
amrex::LPInfo::has_metric_term
bool has_metric_term
Definition AMReX_MLLinOp.H:61
amrex::LinOpEnumType::Location
Location
Definition AMReX_MLLinOp.H:121
amrex::MFInfo
FabArray memory allocation information.
Definition AMReX_FabArray.H:73
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