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
5#include <AMReX_MLCellABecLap.H>
6#include <AMReX_MLALap_K.H>
7#include <AMReX_MultiFabUtil.H>
8
9namespace amrex {
10
18template <typename MF>
20class MLALaplacianT
21 : public MLCellABecLapT<MF>
22{
23public:
24
25 using FAB = typename MF::fab_type;
26 using RT = typename MF::value_type;
27
28 using BCType = LinOpBCType;
29 using Location = typename MLLinOpT<MF>::Location;
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
51 MLALaplacianT (const MLALaplacianT<MF>&) = delete;
52 MLALaplacianT (MLALaplacianT<MF>&&) = delete;
53 MLALaplacianT<MF>& operator= (const MLALaplacianT<MF>&) = delete;
54 MLALaplacianT<MF>& operator= (MLALaplacianT<MF>&&) = delete;
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>
159MLALaplacianT<MF>::MLALaplacianT (const Vector<Geometry>& a_geom,
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
172MLALaplacianT<MF>::define (const Vector<Geometry>& a_geom,
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>
197MLALaplacianT<MF>::~MLALaplacianT () = default;
198
199template <typename MF>
200void
201MLALaplacianT<MF>::setScalars (RT a, RT b) noexcept
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
226MLALaplacianT<MF>::averageDownCoeffs ()
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
243MLALaplacianT<MF>::averageDownCoeffsSameAmrLevel (int amrlev, Vector<MF>& a)
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
264MLALaplacianT<MF>::averageDownCoeffsToCoarseAmrLevel (int flev)
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
279MLALaplacianT<MF>::updateSingularFlag ()
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
338MLALaplacianT<MF>::prepareForSolve ()
339{
340 BL_PROFILE("MLALaplacian::prepareForSolve()");
341 MLCellABecLapT<MF>::prepareForSolve();
342 averageDownCoeffs();
343 updateSingularFlag();
344 m_needs_update = false;
345}
346
347template <typename MF>
348void
349MLALaplacianT<MF>::update ()
350{
351 if (MLCellABecLapT<MF>::needsUpdate()) { MLCellABecLapT<MF>::update(); }
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
367 const GpuArray<RT,AMREX_SPACEDIM> dxinv
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) {
396 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( bx, tbx,
397 {
398 mlalap_adotx_m(tbx, yfab, xfab, afab, dxinv, ascalar, bscalar, dx, probxlo, ncomp);
399 });
400 } else {
401 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( bx, tbx,
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);
412 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( bx2d, tbx2d,
413 {
414 TwoD::mlalap_adotx(tbx2d, yfab2d, xfab2d, afab2d, dhinv, ascalar, bscalar, ncomp);
415 });
416 } else {
417 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( bx, tbx,
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
436 const GpuArray<RT,AMREX_SPACEDIM> dxinv
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) {
464 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( bx, tbx,
465 {
466 mlalap_normalize_m(tbx, fab, afab, dxinv, ascalar, bscalar, dx, probxlo, ncomp);
467 });
468 } else {
469 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( bx, tbx,
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);
479 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( bx2d, tbx2d,
480 {
481 TwoD::mlalap_normalize(tbx2d, fab2d, afab2d, dhinv, ascalar, bscalar, ncomp);
482 });
483 } else {
484 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( bx, tbx,
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 MFItInfo mfi_info;
547 if (Gpu::notInLaunchRegion()) { mfi_info.EnableTiling().SetDynamic(true); }
548
549#ifdef AMREX_USE_OMP
550#pragma omp parallel if (Gpu::notInLaunchRegion())
551#endif
552 for (MFIter mfi(sol,mfi_info); mfi.isValid(); ++mfi)
553 {
554 const auto& m0 = mm0.array(mfi);
555 const auto& m1 = mm1.array(mfi);
556#if (AMREX_SPACEDIM > 1)
557 const auto& m2 = mm2.array(mfi);
558 const auto& m3 = mm3.array(mfi);
559#if (AMREX_SPACEDIM > 2)
560 const auto& m4 = mm4.array(mfi);
561 const auto& m5 = mm5.array(mfi);
562#endif
563#endif
564
565 const Box& tbx = mfi.tilebox();
566 const Box& vbx = mfi.validbox();
567 const auto& solnfab = sol.array(mfi);
568 const auto& rhsfab = rhs.array(mfi);
569 const auto& afab = acoef.array(mfi);
570
571 const auto& f0fab = f0.array(mfi);
572 const auto& f1fab = f1.array(mfi);
573#if (AMREX_SPACEDIM > 1)
574 const auto& f2fab = f2.array(mfi);
575 const auto& f3fab = f3.array(mfi);
576#if (AMREX_SPACEDIM > 2)
577 const auto& f4fab = f4.array(mfi);
578 const auto& f5fab = f5.array(mfi);
579#endif
580#endif
581
582#if (AMREX_SPACEDIM == 1)
583 if (this->m_has_metric_term) {
584 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( tbx, thread_box,
585 {
586 mlalap_gsrb_m(thread_box, solnfab, rhsfab, alpha, dhx,
587 afab,
588 f0fab, m0,
589 f1fab, m1,
590 vbx, redblack,
591 dx, probxlo, ncomp);
592 });
593 } else {
594 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( tbx, thread_box,
595 {
596 mlalap_gsrb(thread_box, solnfab, rhsfab, alpha, dhx,
597 afab,
598 f0fab, m0,
599 f1fab, m1,
600 vbx, redblack, ncomp);
601 });
602 }
603
604#endif
605
606#if (AMREX_SPACEDIM == 2)
607 if (this->m_has_metric_term) {
608 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( tbx, thread_box,
609 {
610 mlalap_gsrb_m(thread_box, solnfab, rhsfab, alpha, dhx, dhy,
611 afab,
612 f0fab, m0,
613 f1fab, m1,
614 f2fab, m2,
615 f3fab, m3,
616 vbx, redblack,
617 dx, probxlo, ncomp);
618 });
619 } else {
620 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( tbx, thread_box,
621 {
622 mlalap_gsrb(thread_box, solnfab, rhsfab, alpha, dhx, dhy,
623 afab,
624 f0fab, m0,
625 f1fab, m1,
626 f2fab, m2,
627 f3fab, m3,
628 vbx, redblack, ncomp);
629 });
630 }
631#endif
632
633#if (AMREX_SPACEDIM == 3)
634 if (this->hasHiddenDimension()) {
635 Box const& tbx_2d = this->compactify(tbx);
636 Box const& vbx_2d = this->compactify(vbx);
637 const auto& solnfab_2d = this->compactify(solnfab);
638 const auto& rhsfab_2d = this->compactify(rhsfab);
639 const auto& afab_2d = this->compactify(afab);
640 const auto& f0fab_2d = this->compactify(this->get_d0(f0fab,f1fab,f2fab));
641 const auto& f1fab_2d = this->compactify(this->get_d1(f0fab,f1fab,f2fab));
642 const auto& f2fab_2d = this->compactify(this->get_d0(f3fab,f4fab,f5fab));
643 const auto& f3fab_2d = this->compactify(this->get_d1(f3fab,f4fab,f5fab));
644 const auto& m0_2d = this->compactify(this->get_d0(m0,m1,m2));
645 const auto& m1_2d = this->compactify(this->get_d1(m0,m1,m2));
646 const auto& m2_2d = this->compactify(this->get_d0(m3,m4,m5));
647 const auto& m3_2d = this->compactify(this->get_d1(m3,m4,m5));
648 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( tbx_2d, thread_box,
649 {
650 TwoD::mlalap_gsrb(thread_box, solnfab_2d, rhsfab_2d, alpha, dh0, dh1,
651 afab_2d,
652 f0fab_2d, m0_2d,
653 f1fab_2d, m1_2d,
654 f2fab_2d, m2_2d,
655 f3fab_2d, m3_2d,
656 vbx_2d, redblack, ncomp);
657 });
658 } else {
659 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( tbx, thread_box,
660 {
661 mlalap_gsrb(thread_box, solnfab, rhsfab, alpha, dhx, dhy, dhz,
662 afab,
663 f0fab, m0,
664 f1fab, m1,
665 f2fab, m2,
666 f3fab, m3,
667 f4fab, m4,
668 f5fab, m5,
669 vbx, redblack, ncomp);
670 });
671 }
672#endif
673 }
674}
675
676template <typename MF>
677void
678MLALaplacianT<MF>::FFlux (int amrlev, const MFIter& mfi,
679 const Array<FAB*,AMREX_SPACEDIM>& flux,
680 const FAB& sol, Location, int face_only) const
681{
682 BL_PROFILE("MLALaplacian::FFlux()");
683
684 const int ncomp = this->getNComp();
685 const int mglev = 0;
686 const Box& box = mfi.tilebox();
687 const Real* dxinv = this->m_geom[amrlev][mglev].InvCellSize();
688
689 AMREX_D_TERM(const auto& fxarr = flux[0]->array();,
690 const auto& fyarr = flux[1]->array();,
691 const auto& fzarr = flux[2]->array(););
692 const auto& solarr = sol.array();
693
694#if (AMREX_SPACEDIM != 3)
695 const RT dx = RT(this->m_geom[amrlev][mglev].CellSize(0));
696 const RT probxlo = RT(this->m_geom[amrlev][mglev].ProbLo(0));
697#endif
698
699#if (AMREX_SPACEDIM == 3)
700 if (face_only) {
701 if (this->hiddenDirection() != 0) {
702 RT fac = m_b_scalar * RT(dxinv[0]);
703 Box blo = amrex::bdryLo(box, 0);
704 int blen = box.length(0);
705 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( blo, tbox,
706 {
707 mlalap_flux_xface(tbox, fxarr, solarr, fac, blen, ncomp);
708 });
709 } else {
710 flux[0]->template setVal<RunOn::Device>(RT(0.0));
711 }
712 if (this->hiddenDirection() != 1) {
713 RT fac = m_b_scalar * RT(dxinv[1]);
714 Box blo = amrex::bdryLo(box, 1);
715 int blen = box.length(1);
716 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( blo, tbox,
717 {
718 mlalap_flux_yface(tbox, fyarr, solarr, fac, blen, ncomp);
719 });
720 } else {
721 flux[1]->template setVal<RunOn::Device>(RT(0.0));
722 }
723 if (this->hiddenDirection() != 2) {
724 RT fac = m_b_scalar * RT(dxinv[2]);
725 Box blo = amrex::bdryLo(box, 2);
726 int blen = box.length(2);
727 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( blo, tbox,
728 {
729 mlalap_flux_zface(tbox, fzarr, solarr, fac, blen, ncomp);
730 });
731 } else {
732 flux[2]->template setVal<RunOn::Device>(RT(0.0));
733 }
734 } else {
735 if (this->hiddenDirection() != 0) {
736 RT fac = m_b_scalar * RT(dxinv[0]);
737 Box bflux = amrex::surroundingNodes(box, 0);
738 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( bflux, tbox,
739 {
740 mlalap_flux_x(tbox, fxarr, solarr, fac, ncomp);
741 });
742 } else {
743 flux[0]->template setVal<RunOn::Device>(RT(0.0));
744 }
745 if (this->hiddenDirection() != 1) {
746 RT fac = m_b_scalar * RT(dxinv[1]);
747 Box bflux = amrex::surroundingNodes(box, 1);
748 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( bflux, tbox,
749 {
750 mlalap_flux_y(tbox, fyarr, solarr, fac, ncomp);
751 });
752 } else {
753 flux[1]->template setVal<RunOn::Device>(RT(0.0));
754 }
755 if (this->hiddenDirection() != 2) {
756 RT fac = m_b_scalar * RT(dxinv[2]);
757 Box bflux = amrex::surroundingNodes(box, 2);
758 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( bflux, tbox,
759 {
760 mlalap_flux_z(tbox, fzarr, solarr, fac, ncomp);
761 });
762 } else {
763 flux[2]->template setVal<RunOn::Device>(RT(0.0));
764 }
765 }
766#elif (AMREX_SPACEDIM == 2)
767 if (face_only) {
768 if (this->hiddenDirection() != 0) {
769 RT fac = m_b_scalar * RT(dxinv[0]);
770 Box blo = amrex::bdryLo(box, 0);
771 int blen = box.length(0);
772 if (this->m_has_metric_term) {
773 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( blo, tbox,
774 {
775 mlalap_flux_xface_m(tbox, fxarr, solarr, fac, blen, dx, probxlo, ncomp);
776 });
777 } else {
778 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( blo, tbox,
779 {
780 mlalap_flux_xface(tbox, fxarr, solarr, fac, blen, ncomp);
781 });
782 }
783 } else {
784 flux[0]->template setVal<RunOn::Device>(RT(0.0));
785 }
786 if (this->hiddenDirection() != 1) {
787 RT fac = m_b_scalar * RT(dxinv[1]);
788 Box blo = amrex::bdryLo(box, 1);
789 int blen = box.length(1);
790 if (this->m_has_metric_term) {
791 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( blo, tbox,
792 {
793 mlalap_flux_yface_m(tbox, fyarr, solarr, fac, blen, dx, probxlo, ncomp);
794 });
795 } else {
796 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( blo, tbox,
797 {
798 mlalap_flux_yface(tbox, fyarr, solarr, fac, blen, ncomp);
799 });
800 }
801 } else {
802 flux[1]->template setVal<RunOn::Device>(RT(0.0));
803 }
804 } else {
805 if (this->hiddenDirection() != 0) {
806 RT fac = m_b_scalar * RT(dxinv[0]);
807 Box bflux = amrex::surroundingNodes(box, 0);
808 if (this->m_has_metric_term) {
809 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( bflux, tbox,
810 {
811 mlalap_flux_x_m(tbox, fxarr, solarr, fac, dx, probxlo, ncomp);
812 });
813 } else {
814 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( bflux, tbox,
815 {
816 mlalap_flux_x(tbox, fxarr, solarr, fac, ncomp);
817 });
818 }
819 } else {
820 flux[0]->template setVal<RunOn::Device>(RT(0.0));
821 }
822 if (this->hiddenDirection() != 1) {
823 RT fac = m_b_scalar * RT(dxinv[1]);
824 Box bflux = amrex::surroundingNodes(box, 1);
825 if (this->m_has_metric_term) {
826 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( bflux, tbox,
827 {
828 mlalap_flux_y_m(tbox, fyarr, solarr, fac, dx, probxlo, ncomp);
829 });
830 } else {
831 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( bflux, tbox,
832 {
833 mlalap_flux_y(tbox, fyarr, solarr, fac, ncomp);
834 });
835 }
836 } else {
837 flux[1]->template setVal<RunOn::Device>(RT(0.0));
838 }
839 }
840#else
841 if (face_only) {
842 RT fac = m_b_scalar * RT(dxinv[0]);
843 Box blo = amrex::bdryLo(box, 0);
844 int blen = box.length(0);
845 if (this->m_has_metric_term) {
846 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( blo, tbox,
847 {
848 mlalap_flux_xface_m(tbox, fxarr, solarr, fac, blen, dx, probxlo, ncomp);
849 });
850 } else {
851 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( blo, tbox,
852 {
853 mlalap_flux_xface(tbox, fxarr, solarr, fac, blen, ncomp);
854 });
855 }
856 } else {
857 RT fac = m_b_scalar * RT(dxinv[0]);
858 Box bflux = amrex::surroundingNodes(box, 0);
859 if (this->m_has_metric_term) {
860 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( bflux, tbox,
861 {
862 mlalap_flux_x_m(tbox, fxarr, solarr, fac, dx, probxlo, ncomp);
863 });
864 } else {
865 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( bflux, tbox,
866 {
867 mlalap_flux_x(tbox, fxarr, solarr, fac, ncomp);
868 });
869 }
870 }
871#endif
872}
873
874extern template class MLALaplacianT<MultiFab>;
875
876using MLALaplacian = MLALaplacianT<MultiFab>;
877
878}
879
880#endif
BL_PROFILE
#define BL_PROFILE(a)
Definition AMReX_BLProfiler.H:562
AMREX_ASSERT
#define AMREX_ASSERT(EX)
Definition AMReX_BLassert.H:38
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_MLCellABecLap.H
AMReX_MultiFabUtil.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::BoxND::length
__host__ __device__ IntVectND< dim > length() const noexcept
Return the length of the BoxND.
Definition AMReX_Box.H:167
amrex::FabFactory
Abstract factory interface for creating, aliasing, and destroying FAB objects.
Definition AMReX_FabFactory.H:73
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::MLALaplacianT::getAScalar
RT getAScalar() const final
Scalar alpha applied to the a term.
Definition AMReX_MLALaplacian.H:112
amrex::MLALaplacianT::operator=
MLALaplacianT< MF > & operator=(const MLALaplacianT< MF > &)=delete
amrex::MLALaplacianT::MLALaplacianT
MLALaplacianT(const MLALaplacianT< MF > &)=delete
amrex::MLALaplacianT::averageDownCoeffsToCoarseAmrLevel
void averageDownCoeffsToCoarseAmrLevel(int flev)
Average a coefficients from fine AMR level flev to flev-1.
Definition AMReX_MLALaplacian.H:264
amrex::MLALaplacianT::~MLALaplacianT
~MLALaplacianT() override
amrex::MLALaplacianT::FAB
typename MF::fab_type FAB
Definition AMReX_MLALaplacian.H:25
amrex::MLALaplacianT::setScalars
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
amrex::MLALaplacianT::FFlux
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:678
amrex::MLALaplacianT::update
void update() override
Update for reuse.
Definition AMReX_MLALaplacian.H:349
amrex::MLALaplacianT::isBottomSingular
bool isBottomSingular() const final
Shortcut for the coarsest level singular flag.
Definition AMReX_MLALaplacian.H:93
amrex::MLALaplacianT::makeNLinOp
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
amrex::MLALaplacianT::getBCoeffs
Array< MF const *, 3 > getBCoeffs(int, int) const final
ALaplacian has no b coefficients; this returns null pointers.
Definition AMReX_MLALaplacian.H:119
amrex::MLALaplacianT::getACoeffs
MF const * getACoeffs(int amrlev, int mglev) const final
Access the stored a coefficient MultiFab for (amrlev,mglev).
Definition AMReX_MLALaplacian.H:116
amrex::MLALaplacianT::prepareForSolve
void prepareForSolve() final
Complete per-level setup (averaging, singularity flags) before solving.
Definition AMReX_MLALaplacian.H:338
amrex::MLALaplacianT::Fapply
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
amrex::MLALaplacianT::averageDownCoeffs
void averageDownCoeffs()
Average a coefficients down across all AMR and MG levels.
Definition AMReX_MLALaplacian.H:226
amrex::MLALaplacianT::MLALaplacianT
MLALaplacianT(MLALaplacianT< MF > &&)=delete
amrex::MLALaplacianT::isSingular
bool isSingular(int amrlev) const final
True if level amrlev is singular.
Definition AMReX_MLALaplacian.H:91
amrex::MLALaplacianT::MLALaplacianT
MLALaplacianT()=default
Construct an empty operator; call define() before use.
amrex::MLALaplacianT::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_MLALaplacian.H:428
amrex::MLALaplacianT::Fsmooth
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
amrex::MLALaplacianT::needsUpdate
bool needsUpdate() const override
Does it need update if it's reused?
Definition AMReX_MLALaplacian.H:83
amrex::MLALaplacianT::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={})
Bind the operator to an AMR hierarchy (no overset support).
Definition AMReX_MLALaplacian.H:172
amrex::MLALaplacianT::getBScalar
RT getBScalar() const final
Scalar beta applied to the Laplacian term.
Definition AMReX_MLALaplacian.H:114
amrex::MLALaplacianT::averageDownCoeffsSameAmrLevel
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
amrex::MLALaplacianT::RT
typename MF::value_type RT
Definition AMReX_MLALaplacian.H:26
amrex::MLALaplacianT::getNComp
int getNComp() const override
Return number of components.
Definition AMReX_MLALaplacian.H:81
amrex::MLALaplacianT::setACoeffs
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
amrex::MLALaplacianT::Location
typename MLLinOpT< MF >::Location Location
Definition AMReX_MLALaplacian.H:29
amrex::MLCellABecLapT
Cell-centered operator that exposes ABec Laplacian helpers to derived classes.
Definition AMReX_MLCellABecLap.H:22
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:309
amrex::MLCellABecLapT::update
void update() override
Average coefficients/metrics when marked dirty.
Definition AMReX_MLCellABecLap.H:302
amrex::MultiMask
Definition AMReX_MultiMask.H:23
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::Vector
This class is a thin wrapper around std::vector. Unlike vector, Vector::operator[] provides bound che...
Definition AMReX_Vector.H:29
amrex::Vector::size
Long size() const noexcept
Definition AMReX_Vector.H:54
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::Array
std::array< T, N > Array
Definition AMReX_Array.H:31
amrex::Gpu::notInLaunchRegion
bool notInLaunchRegion() noexcept
Definition AMReX_GpuControl.H:89
amrex
Definition AMReX_Amr.cpp:50
amrex::average_down
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
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::Abort
void Abort(const std::string &msg)
Print a fatal-error message to stderr and abort execution.
Definition AMReX.cpp:242
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::GpuArray
Fixed-size array that can be used on GPU.
Definition AMReX_Array.H:52
amrex::LPInfo
Configuration knobs for multilevel linear operators (grid agglomeration, metrics, etc....
Definition AMReX_MLLinOp.H:51
amrex::LinOpEnumType::Location
Location
Definition AMReX_MLLinOp.H:119
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