Block-Structured AMR Software Framework
Loading...
Searching...
No Matches
AMReX_MLALaplacian.H
Go to the documentation of this file.
1#ifndef AMREX_MLALAPLACIAN_H_
2#define AMREX_MLALAPLACIAN_H_
3#include <AMReX_Config.H>
4
6#include <AMReX_MLALap_K.H>
8
9namespace amrex {
10
18template <typename MF>
21 : public MLCellABecLapT<MF>
22{
23public:
24
25 using FAB = typename MF::fab_type;
26 using RT = typename MF::value_type;
27
30
32 MLALaplacianT () = default;
43 MLALaplacianT (const Vector<Geometry>& a_geom,
44 const Vector<BoxArray>& a_grids,
45 const Vector<DistributionMapping>& a_dmap,
46 const LPInfo& a_info = LPInfo(),
47 const Vector<FabFactory<FAB> const*>& a_factory = {},
48 int a_ncomp = 1);
49 ~MLALaplacianT () override;
50
55
65 void define (const Vector<Geometry>& a_geom,
66 const Vector<BoxArray>& a_grids,
67 const Vector<DistributionMapping>& a_dmap,
68 const LPInfo& a_info = LPInfo(),
69 const Vector<FabFactory<FAB> const*>& a_factory = {});
70
77 void setScalars (RT a, RT b) noexcept;
79 void setACoeffs (int amrlev, const MF& alpha);
80
81 [[nodiscard]] int getNComp () const override { return m_ncomp; }
82
83 [[nodiscard]] bool needsUpdate () const override {
84 return (m_needs_update || MLCellABecLapT<MF>::needsUpdate());
85 }
86 void update () override;
87
89 void prepareForSolve () final;
91 [[nodiscard]] bool isSingular (int amrlev) const final { return m_is_singular[amrlev]; }
93 [[nodiscard]] bool isBottomSingular () const final { return m_is_singular[0]; }
95 void Fapply (int amrlev, int mglev, MF& out, const MF& in) const final;
97 void Fsmooth (int amrlev, int mglev, MF& sol, const MF& rhs, int redblack) const final;
103 void FFlux (int amrlev, const MFIter& mfi,
104 const Array<FAB*,AMREX_SPACEDIM>& flux,
105 const FAB& sol, Location /* loc */,
106 int face_only=0) const final;
107
109 void normalize (int amrlev, int mglev, MF& mf) const final;
110
112 [[nodiscard]] RT getAScalar () const final { return m_a_scalar; }
114 [[nodiscard]] RT getBScalar () const final { return m_b_scalar; }
116 [[nodiscard]] MF const* getACoeffs (int amrlev, int mglev) const final
117 { return &(m_a_coeffs[amrlev][mglev]); }
119 [[nodiscard]] Array<MF const*,AMREX_SPACEDIM> getBCoeffs (int /*amrlev*/, int /*mglev*/) const final
120 { return {{ AMREX_D_DECL(nullptr,nullptr,nullptr)}}; }
121
122 [[nodiscard]] std::unique_ptr<MLLinOpT<MF>> makeNLinOp (int /*grid_size*/) const final {
123 amrex::Abort("MLALaplacian::makeNLinOp: Not implemented");
124 return std::unique_ptr<MLLinOpT<MF>>{};
125 }
126
133 void averageDownCoeffsSameAmrLevel (int amrlev, Vector<MF>& a);
135 void averageDownCoeffs ();
141 void averageDownCoeffsToCoarseAmrLevel (int flev);
142
143private:
144
145 bool m_needs_update = true;
146
147 RT m_a_scalar = std::numeric_limits<RT>::quiet_NaN();
148 RT m_b_scalar = std::numeric_limits<RT>::quiet_NaN();
149 Vector<Vector<MF> > m_a_coeffs;
150
151 Vector<int> m_is_singular;
152
153 int m_ncomp = 1;
154
155 void updateSingularFlag ();
156};
157
158template <typename MF>
160 const Vector<BoxArray>& a_grids,
161 const Vector<DistributionMapping>& a_dmap,
162 const LPInfo& a_info,
163 const Vector<FabFactory<FAB> const*>& a_factory,
164 int a_ncomp)
165 : m_ncomp(a_ncomp)
166{
167 define(a_geom, a_grids, a_dmap, a_info, a_factory);
168}
169
170template <typename MF>
171void
173 const Vector<BoxArray>& a_grids,
174 const Vector<DistributionMapping>& a_dmap,
175 const LPInfo& a_info,
176 const Vector<FabFactory<FAB> const*>& a_factory)
177{
178 BL_PROFILE("MLALaplacian::define()");
179
180 MLCellABecLapT<MF>::define(a_geom, a_grids, a_dmap, a_info, a_factory);
181
182 const int ncomp = this->getNComp();
183
184 m_a_coeffs.resize(this->m_num_amr_levels);
185 for (int amrlev = 0; amrlev < this->m_num_amr_levels; ++amrlev)
186 {
187 m_a_coeffs[amrlev].resize(this->m_num_mg_levels[amrlev]);
188 for (int mglev = 0; mglev < this->m_num_mg_levels[amrlev]; ++mglev)
189 {
190 m_a_coeffs[amrlev][mglev].define(this->m_grids[amrlev][mglev],
191 this->m_dmap[amrlev][mglev], ncomp, 0);
192 }
193 }
194}
195
196template <typename MF>
198
199template <typename MF>
200void
202{
203 m_a_scalar = a;
204 m_b_scalar = b;
205 if (a == RT(0.0))
206 {
207 for (int amrlev = 0; amrlev < this->m_num_amr_levels; ++amrlev)
208 {
209 m_a_coeffs[amrlev][0].setVal(RT(0.0));
210 }
211 }
212 m_needs_update = true;
213}
214
215template <typename MF>
216void
217MLALaplacianT<MF>::setACoeffs (int amrlev, const MF& alpha)
218{
219 const int ncomp = this->getNComp();
220 m_a_coeffs[amrlev][0].LocalCopy(alpha, 0, 0, ncomp, IntVect(0));
221 m_needs_update = true;
222}
223
224template <typename MF>
225void
227{
228 BL_PROFILE("MLALaplacian::averageDownCoeffs()");
229
230 for (int amrlev = this->m_num_amr_levels-1; amrlev > 0; --amrlev)
231 {
232 auto& fine_a_coeffs = m_a_coeffs[amrlev];
233
234 averageDownCoeffsSameAmrLevel(amrlev, fine_a_coeffs);
235 averageDownCoeffsToCoarseAmrLevel(amrlev);
236 }
237
238 averageDownCoeffsSameAmrLevel(0, m_a_coeffs[0]);
239}
240
241template <typename MF>
242void
244{
245 const int ncomp = this->getNComp();
246 const int nmglevs = a.size();
247 for (int mglev = 1; mglev < nmglevs; ++mglev)
248 {
249 if (m_a_scalar == RT(0.0))
250 {
251 a[mglev].setVal(RT(0.0));
252 }
253 else
254 {
255 AMREX_ASSERT(amrlev == 0 || !this->hasHiddenDimension());
256 IntVect ratio = (amrlev > 0) ? IntVect(this->mg_coarsen_ratio) : this->mg_coarsen_ratio_vec[mglev-1];
257 amrex::average_down(a[mglev-1], a[mglev], 0, ncomp, ratio);
258 }
259 }
260}
261
262template <typename MF>
263void
265{
266 const int ncomp = this->getNComp();
267 auto& fine_a_coeffs = m_a_coeffs[flev ].back();
268 auto& crse_a_coeffs = m_a_coeffs[flev-1].front();
269
270 if (m_a_scalar != RT(0.0)) {
271 // We coarsen from the back of flev to the front of flev-1.
272 // So we use this->mg_coarsen_ratio.
273 amrex::average_down(fine_a_coeffs, crse_a_coeffs, 0, ncomp, this->mg_coarsen_ratio);
274 }
275}
276
277template <typename MF>
278void
280{
281 m_is_singular.clear();
282 m_is_singular.resize(this->m_num_amr_levels, false);
283 auto itlo = std::ranges::find(this->m_lobc[0], BCType::Dirichlet);
284 auto ithi = std::ranges::find(this->m_hibc[0], BCType::Dirichlet);
285 if (itlo == this->m_lobc[0].end() && ithi == this->m_hibc[0].end())
286 { // No Dirichlet
287 for (int alev = 0; alev < this->m_num_amr_levels; ++alev)
288 {
289 if (this->m_domain_covered[alev])
290 {
291 if (m_a_scalar == RT(0.0))
292 {
293 m_is_singular[alev] = true;
294 }
295 else
296 {
297 // We are only testing component 0 here, assuming the others
298 // are similar.
299 RT asum = m_a_coeffs[alev].back().sum(0,IntVect(0));
300 RT amax = m_a_coeffs[alev].back().norminf(0,1,IntVect(0));
301 m_is_singular[alev] = (std::abs(asum) <= amax * RT(1.e-12));
302 }
303 }
304 }
305 }
306
307 if (!m_is_singular[0] && this->m_needs_coarse_data_for_bc &&
308 this->m_coarse_fine_bc_type == BCType::Neumann)
309 {
310 bool lev0_a_is_zero = false;
311 if (m_a_scalar == RT(0.0)) {
312 lev0_a_is_zero = true;
313 } else {
314 RT asum = m_a_coeffs[0].back().sum(0,IntVect(0));
315 RT amax = m_a_coeffs[0].back().norminf(0,1,IntVect(0));
316 lev0_a_is_zero = std::abs(asum) <= amax * RT(1.e-12);
317 }
318
319 if (lev0_a_is_zero) {
320 auto bbox = this->m_grids[0][0].minimalBox();
321 for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
322 if (this->m_lobc[0][idim] == BCType::Dirichlet) {
323 bbox.growLo(idim,1);
324 }
325 if (this->m_hibc[0][idim] == BCType::Dirichlet) {
326 bbox.growHi(idim,1);
327 }
328 }
329 if (this->m_geom[0][0].Domain().contains(bbox)) {
330 m_is_singular[0] = true;
331 }
332 }
333 }
334}
335
336template <typename MF>
337void
339{
340 BL_PROFILE("MLALaplacian::prepareForSolve()");
342 averageDownCoeffs();
343 updateSingularFlag();
344 m_needs_update = false;
345}
346
347template <typename MF>
348void
350{
352 averageDownCoeffs();
353 updateSingularFlag();
354 m_needs_update = false;
355}
356
357template <typename MF>
358void
359MLALaplacianT<MF>::Fapply (int amrlev, int mglev, MF& out, const MF& in) const
360{
361 BL_PROFILE("MLALaplacian::Fapply()");
362
363 const int ncomp = this->getNComp();
364
365 const MF& acoef = m_a_coeffs[amrlev][mglev];
366
368 {AMREX_D_DECL(RT(this->m_geom[amrlev][mglev].InvCellSize(0)),
369 RT(this->m_geom[amrlev][mglev].InvCellSize(1)),
370 RT(this->m_geom[amrlev][mglev].InvCellSize(2)))};
371#if (AMREX_SPACEDIM < 3)
372 const RT dx = RT(this->m_geom[amrlev][mglev].CellSize(0));
373 const RT probxlo = RT(this->m_geom[amrlev][mglev].ProbLo(0));
374#endif
375
376#if (AMREX_SPACEDIM == 3)
377 GpuArray<RT,2> dhinv {this->get_d0(dxinv[0], dxinv[1], dxinv[2]),
378 this->get_d1(dxinv[0], dxinv[1], dxinv[2])};
379#endif
380
381 const RT ascalar = m_a_scalar;
382 const RT bscalar = m_b_scalar;
383
384#ifdef AMREX_USE_OMP
385#pragma omp parallel if (Gpu::notInLaunchRegion())
386#endif
387 for (MFIter mfi(out, TilingIfNotGPU()); mfi.isValid(); ++mfi)
388 {
389 const Box& bx = mfi.tilebox();
390 const auto& xfab = in.array(mfi);
391 const auto& yfab = out.array(mfi);
392 const auto& afab = acoef.array(mfi);
393
394#if (AMREX_SPACEDIM != 3)
395 if (this->m_has_metric_term) {
397 {
398 mlalap_adotx_m(tbx, yfab, xfab, afab, dxinv, ascalar, bscalar, dx, probxlo, ncomp);
399 });
400 } else {
402 {
403 mlalap_adotx(tbx, yfab, xfab, afab, dxinv, ascalar, bscalar, ncomp);
404 });
405 }
406#else
407 if (this->hasHiddenDimension()) {
408 Box const& bx2d = this->compactify(bx);
409 const auto& xfab2d = this->compactify(xfab);
410 const auto& yfab2d = this->compactify(yfab);
411 const auto& afab2d = this->compactify(afab);
413 {
414 TwoD::mlalap_adotx(tbx2d, yfab2d, xfab2d, afab2d, dhinv, ascalar, bscalar, ncomp);
415 });
416 } else {
418 {
419 mlalap_adotx(tbx, yfab, xfab, afab, dxinv, ascalar, bscalar, ncomp);
420 });
421 }
422#endif
423 }
424}
425
426template <typename MF>
427void
428MLALaplacianT<MF>::normalize (int amrlev, int mglev, MF& mf) const
429{
430 BL_PROFILE("MLALaplacian::normalize()");
431
432 const int ncomp = this->getNComp();
433
434 const MF& acoef = m_a_coeffs[amrlev][mglev];
435
437 {AMREX_D_DECL(RT(this->m_geom[amrlev][mglev].InvCellSize(0)),
438 RT(this->m_geom[amrlev][mglev].InvCellSize(1)),
439 RT(this->m_geom[amrlev][mglev].InvCellSize(2)))};
440#if (AMREX_SPACEDIM < 3)
441 const RT dx = RT(this->m_geom[amrlev][mglev].CellSize(0));
442 const RT probxlo = RT(this->m_geom[amrlev][mglev].ProbLo(0));
443#endif
444
445#if (AMREX_SPACEDIM == 3)
446 GpuArray<RT,2> dhinv {this->get_d0(dxinv[0], dxinv[1], dxinv[2]),
447 this->get_d1(dxinv[0], dxinv[1], dxinv[2])};
448#endif
449
450 const RT ascalar = m_a_scalar;
451 const RT bscalar = m_b_scalar;
452
453#ifdef AMREX_USE_OMP
454#pragma omp parallel if (Gpu::notInLaunchRegion())
455#endif
456 for (MFIter mfi(mf, TilingIfNotGPU()); mfi.isValid(); ++mfi)
457 {
458 const Box& bx = mfi.tilebox();
459 const auto& fab = mf.array(mfi);
460 const auto& afab = acoef.array(mfi);
461
462#if (AMREX_SPACEDIM != 3)
463 if (this->m_has_metric_term) {
465 {
466 mlalap_normalize_m(tbx, fab, afab, dxinv, ascalar, bscalar, dx, probxlo, ncomp);
467 });
468 } else {
470 {
471 mlalap_normalize(tbx, fab, afab, dxinv, ascalar, bscalar, ncomp);
472 });
473 }
474#else
475 if (this->hasHiddenDimension()) {
476 Box const& bx2d = this->compactify(bx);
477 const auto& fab2d = this->compactify(fab);
478 const auto& afab2d = this->compactify(afab);
480 {
481 TwoD::mlalap_normalize(tbx2d, fab2d, afab2d, dhinv, ascalar, bscalar, ncomp);
482 });
483 } else {
485 {
486 mlalap_normalize(tbx, fab, afab, dxinv, ascalar, bscalar, ncomp);
487 });
488 }
489#endif
490 }
491}
492
493template <typename MF>
494void
495MLALaplacianT<MF>::Fsmooth (int amrlev, int mglev, MF& sol, const MF& rhs, int redblack) const
496{
497 BL_PROFILE("MLALaplacian::Fsmooth()");
498
499 const int ncomp = this->getNComp();
500
501 const MF& acoef = m_a_coeffs[amrlev][mglev];
502 const auto& undrrelxr = this->m_undrrelxr[amrlev][mglev];
503 const auto& maskvals = this->m_maskvals [amrlev][mglev];
504
505 OrientationIter oitr;
506
507 const auto& f0 = undrrelxr[oitr()]; ++oitr;
508 const auto& f1 = undrrelxr[oitr()]; ++oitr;
509#if (AMREX_SPACEDIM > 1)
510 const auto& f2 = undrrelxr[oitr()]; ++oitr;
511 const auto& f3 = undrrelxr[oitr()]; ++oitr;
512#if (AMREX_SPACEDIM > 2)
513 const auto& f4 = undrrelxr[oitr()]; ++oitr;
514 const auto& f5 = undrrelxr[oitr()]; ++oitr;
515#endif
516#endif
517
518 const MultiMask& mm0 = maskvals[0];
519 const MultiMask& mm1 = maskvals[1];
520#if (AMREX_SPACEDIM > 1)
521 const MultiMask& mm2 = maskvals[2];
522 const MultiMask& mm3 = maskvals[3];
523#if (AMREX_SPACEDIM > 2)
524 const MultiMask& mm4 = maskvals[4];
525 const MultiMask& mm5 = maskvals[5];
526#endif
527#endif
528
529 const Real* dxinv = this->m_geom[amrlev][mglev].InvCellSize();
530 AMREX_D_TERM(const RT dhx = m_b_scalar*RT(dxinv[0]*dxinv[0]);,
531 const RT dhy = m_b_scalar*RT(dxinv[1]*dxinv[1]);,
532 const RT dhz = m_b_scalar*RT(dxinv[2]*dxinv[2]););
533
534#if (AMREX_SPACEDIM == 3)
535 RT dh0 = this->get_d0(dhx, dhy, dhz);
536 RT dh1 = this->get_d1(dhx, dhy, dhz);
537#endif
538
539#if (AMREX_SPACEDIM < 3)
540 const RT dx = RT(this->m_geom[amrlev][mglev].CellSize(0));
541 const RT probxlo = RT(this->m_geom[amrlev][mglev].ProbLo(0));
542#endif
543
544 const RT alpha = m_a_scalar;
545
546#ifdef AMREX_USE_GPU
547 if (Gpu::inLaunchRegion() && !this->hasHiddenDimension())
548 {
549 const auto& m0ma = mm0.const_arrays();
550 const auto& m1ma = mm1.const_arrays();
551 const auto& f0ma = f0.const_arrays();
552 const auto& f1ma = f1.const_arrays();
553#if (AMREX_SPACEDIM > 1)
554 const auto& m2ma = mm2.const_arrays();
555 const auto& m3ma = mm3.const_arrays();
556 const auto& f2ma = f2.const_arrays();
557 const auto& f3ma = f3.const_arrays();
558#if (AMREX_SPACEDIM > 2)
559 const auto& m4ma = mm4.const_arrays();
560 const auto& m5ma = mm5.const_arrays();
561 const auto& f4ma = f4.const_arrays();
562 const auto& f5ma = f5.const_arrays();
563#endif
564#endif
565 const auto& solnma = sol.arrays();
566 const auto& rhsma = rhs.const_arrays();
567 const auto& ama = acoef.const_arrays();
568
569#if (AMREX_SPACEDIM < 3)
570 if (this->m_has_metric_term) {
571 ParallelForRedBlack(sol, ncomp, redblack,
572 [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k, int n) noexcept
573 {
574 Box vbx(ama[box_no]);
575 mlalap_gsrb_m(i,j,k,n, solnma[box_no], rhsma[box_no], alpha,
576 AMREX_D_DECL(dhx, dhy, dhz), ama[box_no],
577 f0ma[box_no], m0ma[box_no],
578 f1ma[box_no], m1ma[box_no],
579#if (AMREX_SPACEDIM > 1)
580 f2ma[box_no], m2ma[box_no],
581 f3ma[box_no], m3ma[box_no],
582#endif
583 vbx, redblack, dx, probxlo);
584 });
585 } else
586#endif
587 {
588 ParallelForRedBlack(sol, ncomp, redblack,
589 [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k, int n) noexcept
590 {
591 Box vbx(ama[box_no]);
592 mlalap_gsrb(i,j,k,n, solnma[box_no], rhsma[box_no], alpha,
593 AMREX_D_DECL(dhx, dhy, dhz), ama[box_no],
594 f0ma[box_no], m0ma[box_no],
595 f1ma[box_no], m1ma[box_no],
596#if (AMREX_SPACEDIM > 1)
597 f2ma[box_no], m2ma[box_no],
598 f3ma[box_no], m3ma[box_no],
599#if (AMREX_SPACEDIM > 2)
600 f4ma[box_no], m4ma[box_no],
601 f5ma[box_no], m5ma[box_no],
602#endif
603#endif
604 vbx, redblack);
605 });
606 }
607 if (!Gpu::inNoSyncRegion()) {
609 }
610 return;
611 }
612#endif
613
614 MFItInfo mfi_info;
615 if (Gpu::notInLaunchRegion()) { mfi_info.EnableTiling().SetDynamic(true); }
616
617#ifdef AMREX_USE_OMP
618#pragma omp parallel if (Gpu::notInLaunchRegion())
619#endif
620 for (MFIter mfi(sol,mfi_info); mfi.isValid(); ++mfi)
621 {
622 const auto& m0 = mm0.array(mfi);
623 const auto& m1 = mm1.array(mfi);
624#if (AMREX_SPACEDIM > 1)
625 const auto& m2 = mm2.array(mfi);
626 const auto& m3 = mm3.array(mfi);
627#if (AMREX_SPACEDIM > 2)
628 const auto& m4 = mm4.array(mfi);
629 const auto& m5 = mm5.array(mfi);
630#endif
631#endif
632
633 const Box& tbx = mfi.tilebox();
634 const Box& vbx = mfi.validbox();
635 const auto& solnfab = sol.array(mfi);
636 const auto& rhsfab = rhs.array(mfi);
637 const auto& afab = acoef.array(mfi);
638
639 const auto& f0fab = f0.array(mfi);
640 const auto& f1fab = f1.array(mfi);
641#if (AMREX_SPACEDIM > 1)
642 const auto& f2fab = f2.array(mfi);
643 const auto& f3fab = f3.array(mfi);
644#if (AMREX_SPACEDIM > 2)
645 const auto& f4fab = f4.array(mfi);
646 const auto& f5fab = f5.array(mfi);
647#endif
648#endif
649
650#if (AMREX_SPACEDIM == 1)
651 if (this->m_has_metric_term) {
653 {
654 mlalap_gsrb_m(thread_box, solnfab, rhsfab, alpha, dhx,
655 afab,
656 f0fab, m0,
657 f1fab, m1,
658 vbx, redblack,
659 dx, probxlo, ncomp);
660 });
661 } else {
663 {
664 mlalap_gsrb(thread_box, solnfab, rhsfab, alpha, dhx,
665 afab,
666 f0fab, m0,
667 f1fab, m1,
668 vbx, redblack, ncomp);
669 });
670 }
671
672#endif
673
674#if (AMREX_SPACEDIM == 2)
675 if (this->m_has_metric_term) {
677 {
678 mlalap_gsrb_m(thread_box, solnfab, rhsfab, alpha, dhx, dhy,
679 afab,
680 f0fab, m0,
681 f1fab, m1,
682 f2fab, m2,
683 f3fab, m3,
684 vbx, redblack,
685 dx, probxlo, ncomp);
686 });
687 } else {
689 {
690 mlalap_gsrb(thread_box, solnfab, rhsfab, alpha, dhx, dhy,
691 afab,
692 f0fab, m0,
693 f1fab, m1,
694 f2fab, m2,
695 f3fab, m3,
696 vbx, redblack, ncomp);
697 });
698 }
699#endif
700
701#if (AMREX_SPACEDIM == 3)
702 if (this->hasHiddenDimension()) {
703 Box const& tbx_2d = this->compactify(tbx);
704 Box const& vbx_2d = this->compactify(vbx);
705 const auto& solnfab_2d = this->compactify(solnfab);
706 const auto& rhsfab_2d = this->compactify(rhsfab);
707 const auto& afab_2d = this->compactify(afab);
708 const auto& f0fab_2d = this->compactify(this->get_d0(f0fab,f1fab,f2fab));
709 const auto& f1fab_2d = this->compactify(this->get_d1(f0fab,f1fab,f2fab));
710 const auto& f2fab_2d = this->compactify(this->get_d0(f3fab,f4fab,f5fab));
711 const auto& f3fab_2d = this->compactify(this->get_d1(f3fab,f4fab,f5fab));
712 const auto& m0_2d = this->compactify(this->get_d0(m0,m1,m2));
713 const auto& m1_2d = this->compactify(this->get_d1(m0,m1,m2));
714 const auto& m2_2d = this->compactify(this->get_d0(m3,m4,m5));
715 const auto& m3_2d = this->compactify(this->get_d1(m3,m4,m5));
717 {
718 TwoD::mlalap_gsrb(thread_box, solnfab_2d, rhsfab_2d, alpha, dh0, dh1,
719 afab_2d,
720 f0fab_2d, m0_2d,
721 f1fab_2d, m1_2d,
722 f2fab_2d, m2_2d,
723 f3fab_2d, m3_2d,
724 vbx_2d, redblack, ncomp);
725 });
726 } else {
728 {
729 mlalap_gsrb(thread_box, solnfab, rhsfab, alpha, dhx, dhy, dhz,
730 afab,
731 f0fab, m0,
732 f1fab, m1,
733 f2fab, m2,
734 f3fab, m3,
735 f4fab, m4,
736 f5fab, m5,
737 vbx, redblack, ncomp);
738 });
739 }
740#endif
741 }
742}
743
744template <typename MF>
745void
746MLALaplacianT<MF>::FFlux (int amrlev, const MFIter& mfi,
747 const Array<FAB*,AMREX_SPACEDIM>& flux,
748 const FAB& sol, Location, int face_only) const
749{
750 BL_PROFILE("MLALaplacian::FFlux()");
751
752 const int ncomp = this->getNComp();
753 const int mglev = 0;
754 const Box& box = mfi.tilebox();
755 const Real* dxinv = this->m_geom[amrlev][mglev].InvCellSize();
756
757 AMREX_D_TERM(const auto& fxarr = flux[0]->array();,
758 const auto& fyarr = flux[1]->array();,
759 const auto& fzarr = flux[2]->array(););
760 const auto& solarr = sol.array();
761
762#if (AMREX_SPACEDIM != 3)
763 const RT dx = RT(this->m_geom[amrlev][mglev].CellSize(0));
764 const RT probxlo = RT(this->m_geom[amrlev][mglev].ProbLo(0));
765#endif
766
767#if (AMREX_SPACEDIM == 3)
768 if (face_only) {
769 if (this->hiddenDirection() != 0) {
770 RT fac = m_b_scalar * RT(dxinv[0]);
771 Box blo = amrex::bdryLo(box, 0);
772 int blen = box.length(0);
774 {
775 mlalap_flux_xface(tbox, fxarr, solarr, fac, blen, ncomp);
776 });
777 } else {
778 flux[0]->template setVal<RunOn::Device>(RT(0.0));
779 }
780 if (this->hiddenDirection() != 1) {
781 RT fac = m_b_scalar * RT(dxinv[1]);
782 Box blo = amrex::bdryLo(box, 1);
783 int blen = box.length(1);
785 {
786 mlalap_flux_yface(tbox, fyarr, solarr, fac, blen, ncomp);
787 });
788 } else {
789 flux[1]->template setVal<RunOn::Device>(RT(0.0));
790 }
791 if (this->hiddenDirection() != 2) {
792 RT fac = m_b_scalar * RT(dxinv[2]);
793 Box blo = amrex::bdryLo(box, 2);
794 int blen = box.length(2);
796 {
797 mlalap_flux_zface(tbox, fzarr, solarr, fac, blen, ncomp);
798 });
799 } else {
800 flux[2]->template setVal<RunOn::Device>(RT(0.0));
801 }
802 } else {
803 if (this->hiddenDirection() != 0) {
804 RT fac = m_b_scalar * RT(dxinv[0]);
805 Box bflux = amrex::surroundingNodes(box, 0);
807 {
808 mlalap_flux_x(tbox, fxarr, solarr, fac, ncomp);
809 });
810 } else {
811 flux[0]->template setVal<RunOn::Device>(RT(0.0));
812 }
813 if (this->hiddenDirection() != 1) {
814 RT fac = m_b_scalar * RT(dxinv[1]);
815 Box bflux = amrex::surroundingNodes(box, 1);
817 {
818 mlalap_flux_y(tbox, fyarr, solarr, fac, ncomp);
819 });
820 } else {
821 flux[1]->template setVal<RunOn::Device>(RT(0.0));
822 }
823 if (this->hiddenDirection() != 2) {
824 RT fac = m_b_scalar * RT(dxinv[2]);
825 Box bflux = amrex::surroundingNodes(box, 2);
827 {
828 mlalap_flux_z(tbox, fzarr, solarr, fac, ncomp);
829 });
830 } else {
831 flux[2]->template setVal<RunOn::Device>(RT(0.0));
832 }
833 }
834#elif (AMREX_SPACEDIM == 2)
835 if (face_only) {
836 if (this->hiddenDirection() != 0) {
837 RT fac = m_b_scalar * RT(dxinv[0]);
838 Box blo = amrex::bdryLo(box, 0);
839 int blen = box.length(0);
840 if (this->m_has_metric_term) {
842 {
843 mlalap_flux_xface_m(tbox, fxarr, solarr, fac, blen, dx, probxlo, ncomp);
844 });
845 } else {
847 {
848 mlalap_flux_xface(tbox, fxarr, solarr, fac, blen, ncomp);
849 });
850 }
851 } else {
852 flux[0]->template setVal<RunOn::Device>(RT(0.0));
853 }
854 if (this->hiddenDirection() != 1) {
855 RT fac = m_b_scalar * RT(dxinv[1]);
856 Box blo = amrex::bdryLo(box, 1);
857 int blen = box.length(1);
858 if (this->m_has_metric_term) {
860 {
861 mlalap_flux_yface_m(tbox, fyarr, solarr, fac, blen, dx, probxlo, ncomp);
862 });
863 } else {
865 {
866 mlalap_flux_yface(tbox, fyarr, solarr, fac, blen, ncomp);
867 });
868 }
869 } else {
870 flux[1]->template setVal<RunOn::Device>(RT(0.0));
871 }
872 } else {
873 if (this->hiddenDirection() != 0) {
874 RT fac = m_b_scalar * RT(dxinv[0]);
875 Box bflux = amrex::surroundingNodes(box, 0);
876 if (this->m_has_metric_term) {
878 {
879 mlalap_flux_x_m(tbox, fxarr, solarr, fac, dx, probxlo, ncomp);
880 });
881 } else {
883 {
884 mlalap_flux_x(tbox, fxarr, solarr, fac, ncomp);
885 });
886 }
887 } else {
888 flux[0]->template setVal<RunOn::Device>(RT(0.0));
889 }
890 if (this->hiddenDirection() != 1) {
891 RT fac = m_b_scalar * RT(dxinv[1]);
892 Box bflux = amrex::surroundingNodes(box, 1);
893 if (this->m_has_metric_term) {
895 {
896 mlalap_flux_y_m(tbox, fyarr, solarr, fac, dx, probxlo, ncomp);
897 });
898 } else {
900 {
901 mlalap_flux_y(tbox, fyarr, solarr, fac, ncomp);
902 });
903 }
904 } else {
905 flux[1]->template setVal<RunOn::Device>(RT(0.0));
906 }
907 }
908#else
909 if (face_only) {
910 RT fac = m_b_scalar * RT(dxinv[0]);
911 Box blo = amrex::bdryLo(box, 0);
912 int blen = box.length(0);
913 if (this->m_has_metric_term) {
915 {
916 mlalap_flux_xface_m(tbox, fxarr, solarr, fac, blen, dx, probxlo, ncomp);
917 });
918 } else {
920 {
921 mlalap_flux_xface(tbox, fxarr, solarr, fac, blen, ncomp);
922 });
923 }
924 } else {
925 RT fac = m_b_scalar * RT(dxinv[0]);
926 Box bflux = amrex::surroundingNodes(box, 0);
927 if (this->m_has_metric_term) {
929 {
930 mlalap_flux_x_m(tbox, fxarr, solarr, fac, dx, probxlo, ncomp);
931 });
932 } else {
934 {
935 mlalap_flux_x(tbox, fxarr, solarr, fac, ncomp);
936 });
937 }
938 }
939#endif
940}
941
942extern template class MLALaplacianT<MultiFab>;
943
945
946}
947
948#endif
#define BL_PROFILE(a)
Definition AMReX_BLProfiler.H:562
#define AMREX_ASSERT(EX)
Definition AMReX_BLassert.H:38
#define AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE(TN, TI, block)
Definition AMReX_GpuLaunchMacrosC.nolint.H:4
#define AMREX_GPU_DEVICE
Definition AMReX_GpuQualifiers.H:18
#define AMREX_D_TERM(a, b, c)
Definition AMReX_SPACE.H:172
#define AMREX_D_DECL(a, b, c)
Definition AMReX_SPACE.H:171
__host__ __device__ IntVectND< dim > length() const noexcept
Return the length of the BoxND.
Definition AMReX_Box.H:167
Abstract factory interface for creating, aliasing, and destroying FAB objects.
Definition AMReX_FabFactory.H:73
Iterator for looping ever tiles and boxes of amrex::FabArray based containers.
Definition AMReX_MFIter.H:88
Box tilebox() const noexcept
Return the tile Box at the current index.
Definition AMReX_MFIter.cpp:389
bool isValid() const noexcept
Is the iterator valid i.e. is it associated with a FAB?
Definition AMReX_MFIter.H:176
Multi-component ALaplacian (a scalar plus optional spatial a coeffs).
Definition AMReX_MLALaplacian.H:22
RT getAScalar() const final
Scalar alpha applied to the a term.
Definition AMReX_MLALaplacian.H:112
MLALaplacianT< MF > & operator=(const MLALaplacianT< MF > &)=delete
MLALaplacianT(const MLALaplacianT< MF > &)=delete
void averageDownCoeffsToCoarseAmrLevel(int flev)
Average a coefficients from fine AMR level flev to flev-1.
Definition AMReX_MLALaplacian.H:264
~MLALaplacianT() override
typename MF::fab_type FAB
Definition AMReX_MLALaplacian.H:25
void setScalars(RT a, RT b) noexcept
Set constant scalars a and b in a \phi - b \nabla^2 \phi.
Definition AMReX_MLALaplacian.H:201
void FFlux(int amrlev, const MFIter &mfi, const Array< FAB *, 3 > &flux, const FAB &sol, Location, int face_only=0) const final
Produce face fluxes on AMR level amrlev for the tilebox described by mfi using sol,...
Definition AMReX_MLALaplacian.H:746
void update() override
Update for reuse.
Definition AMReX_MLALaplacian.H:349
bool isBottomSingular() const final
Shortcut for the coarsest level singular flag.
Definition AMReX_MLALaplacian.H:93
std::unique_ptr< MLLinOpT< MF > > makeNLinOp(int) const final
Create the NSolve counterpart of this operator with the requested grid size.
Definition AMReX_MLALaplacian.H:122
Array< MF const *, 3 > getBCoeffs(int, int) const final
ALaplacian has no b coefficients; this returns null pointers.
Definition AMReX_MLALaplacian.H:119
MF const * getACoeffs(int amrlev, int mglev) const final
Access the stored a coefficient MultiFab for (amrlev,mglev).
Definition AMReX_MLALaplacian.H:116
void prepareForSolve() final
Complete per-level setup (averaging, singularity flags) before solving.
Definition AMReX_MLALaplacian.H:338
void Fapply(int amrlev, int mglev, MF &out, const MF &in) const final
Apply the ALaplacian to in (writing out) on (amrlev,mglev).
Definition AMReX_MLALaplacian.H:359
void averageDownCoeffs()
Average a coefficients down across all AMR and MG levels.
Definition AMReX_MLALaplacian.H:226
MLALaplacianT(MLALaplacianT< MF > &&)=delete
bool isSingular(int amrlev) const final
True if level amrlev is singular.
Definition AMReX_MLALaplacian.H:91
MLALaplacianT()=default
Construct an empty operator; call define() before use.
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_MLALaplacian.H:428
void Fsmooth(int amrlev, int mglev, MF &sol, const MF &rhs, int redblack) const final
Run a smoothing sweep on (amrlev,mglev). redblack selects the red (0) or black (1) half of the grid.
Definition AMReX_MLALaplacian.H:495
bool needsUpdate() const override
Does it need update if it's reused?
Definition AMReX_MLALaplacian.H:83
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={})
Bind the operator to an AMR hierarchy (no overset support).
Definition AMReX_MLALaplacian.H:172
RT getBScalar() const final
Scalar beta applied to the Laplacian term.
Definition AMReX_MLALaplacian.H:114
void averageDownCoeffsSameAmrLevel(int amrlev, Vector< MF > &a)
Average a coefficients down within a single AMR level (fine-to-coarse MG).
Definition AMReX_MLALaplacian.H:243
typename MF::value_type RT
Definition AMReX_MLALaplacian.H:26
int getNComp() const override
Return number of components.
Definition AMReX_MLALaplacian.H:81
void setACoeffs(int amrlev, const MF &alpha)
Provide per-cell a coefficients on AMR level amrlev (stored directly in alpha).
Definition AMReX_MLALaplacian.H:217
typename MLLinOpT< MF >::Location Location
Definition AMReX_MLALaplacian.H:29
Cell-centered operator that exposes ABec Laplacian helpers to derived classes.
Definition AMReX_MLCellABecLap.H:22
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
void prepareForSolve() override
Standard hook called before MLMG iterates (fixes BC data, etc.).
Definition AMReX_MLCellABecLap.H:319
void update() override
Average coefficients/metrics when marked dirty.
Definition AMReX_MLCellABecLap.H:312
Definition AMReX_MultiMask.H:23
MultiArray4< int const > const_arrays() const noexcept
Return const multi-array views (alias of arrays()).
Definition AMReX_MultiMask.H:86
Array4< int const > array(const MFIter &mfi) const noexcept
Return an Array4 view (const) for iterator mfi.
Definition AMReX_MultiMask.H:69
An Iterator over the Orientation of Faces of a Box.
Definition AMReX_Orientation.H:135
This class is a thin wrapper around std::vector. Unlike vector, Vector::operator[] provides bound che...
Definition AMReX_Vector.H:29
Long size() const noexcept
Definition AMReX_Vector.H:54
amrex_real Real
Floating Point Type for Fields.
Definition AMReX_REAL.H:80
__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
__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
std::array< T, N > Array
Definition AMReX_Array.H:31
void streamSynchronize() noexcept
Definition AMReX_GpuDevice.H:310
bool inLaunchRegion() noexcept
Definition AMReX_GpuControl.H:88
bool notInLaunchRegion() noexcept
Definition AMReX_GpuControl.H:89
bool inNoSyncRegion() noexcept
Definition AMReX_GpuControl.H:148
Definition AMReX_Amr.cpp:50
void average_down(const MultiFab &S_fine, MultiFab &S_crse, const Geometry &fgeom, const Geometry &cgeom, int scomp, int ncomp, int rr)
Definition AMReX_MultiFabUtil.cpp:359
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
LinOpBCType
Definition AMReX_LO_BCTYPES.H:27
IntVectND< 3 > IntVect
IntVect is an alias for amrex::IntVectND instantiated with AMREX_SPACEDIM.
Definition AMReX_BaseFwd.H:38
bool TilingIfNotGPU() noexcept
Definition AMReX_MFIter.H:12
void Abort(const std::string &msg)
Print a fatal-error message to stderr and abort execution.
Definition AMReX.cpp:244
__host__ __device__ Dim3 end(BoxND< dim > const &box) noexcept
Return the iterator end coordinate of box as Dim3.
Definition AMReX_Box.H:2257
Fixed-size array that can be used on GPU.
Definition AMReX_Array.H:52
Configuration knobs for multilevel linear operators (grid agglomeration, metrics, etc....
Definition AMReX_MLLinOp.H:53
Location
Definition AMReX_MLLinOp.H:121
Definition AMReX_MFIter.H:20
MFItInfo & SetDynamic(bool f) noexcept
Definition AMReX_MFIter.H:43
MFItInfo & EnableTiling(const IntVect &ts=FabArrayBase::mfiter_tile_size) noexcept
Definition AMReX_MFIter.H:31