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}
213
214template <typename MF>
215void
216MLALaplacianT<MF>::setACoeffs (int amrlev, const MF& alpha)
217{
218 const int ncomp = this->getNComp();
219 m_a_coeffs[amrlev][0].LocalCopy(alpha, 0, 0, ncomp, IntVect(0));
220 m_needs_update = true;
221}
222
223template <typename MF>
224void
225MLALaplacianT<MF>::averageDownCoeffs ()
226{
227 BL_PROFILE("MLALaplacian::averageDownCoeffs()");
228
229 for (int amrlev = this->m_num_amr_levels-1; amrlev > 0; --amrlev)
230 {
231 auto& fine_a_coeffs = m_a_coeffs[amrlev];
232
233 averageDownCoeffsSameAmrLevel(amrlev, fine_a_coeffs);
234 averageDownCoeffsToCoarseAmrLevel(amrlev);
235 }
236
237 averageDownCoeffsSameAmrLevel(0, m_a_coeffs[0]);
238}
239
240template <typename MF>
241void
242MLALaplacianT<MF>::averageDownCoeffsSameAmrLevel (int amrlev, Vector<MF>& a)
243{
244 const int ncomp = this->getNComp();
245 const int nmglevs = a.size();
246 for (int mglev = 1; mglev < nmglevs; ++mglev)
247 {
248 if (m_a_scalar == RT(0.0))
249 {
250 a[mglev].setVal(RT(0.0));
251 }
252 else
253 {
254 AMREX_ASSERT(amrlev == 0 || !this->hasHiddenDimension());
255 IntVect ratio = (amrlev > 0) ? IntVect(this->mg_coarsen_ratio) : this->mg_coarsen_ratio_vec[mglev-1];
256 amrex::average_down(a[mglev-1], a[mglev], 0, ncomp, ratio);
257 }
258 }
259}
260
261template <typename MF>
262void
263MLALaplacianT<MF>::averageDownCoeffsToCoarseAmrLevel (int flev)
264{
265 const int ncomp = this->getNComp();
266 auto& fine_a_coeffs = m_a_coeffs[flev ].back();
267 auto& crse_a_coeffs = m_a_coeffs[flev-1].front();
268
269 if (m_a_scalar != RT(0.0)) {
270 // We coarsen from the back of flev to the front of flev-1.
271 // So we use this->mg_coarsen_ratio.
272 amrex::average_down(fine_a_coeffs, crse_a_coeffs, 0, ncomp, this->mg_coarsen_ratio);
273 }
274}
275
276template <typename MF>
277void
278MLALaplacianT<MF>::updateSingularFlag ()
279{
280 m_is_singular.clear();
281 m_is_singular.resize(this->m_num_amr_levels, false);
282 auto itlo = std::ranges::find(this->m_lobc[0], BCType::Dirichlet);
283 auto ithi = std::ranges::find(this->m_hibc[0], BCType::Dirichlet);
284 if (itlo == this->m_lobc[0].end() && ithi == this->m_hibc[0].end())
285 { // No Dirichlet
286 for (int alev = 0; alev < this->m_num_amr_levels; ++alev)
287 {
288 if (this->m_domain_covered[alev])
289 {
290 if (m_a_scalar == RT(0.0))
291 {
292 m_is_singular[alev] = true;
293 }
294 else
295 {
296 // We are only testing component 0 here, assuming the others
297 // are similar.
298 RT asum = m_a_coeffs[alev].back().sum(0,IntVect(0));
299 RT amax = m_a_coeffs[alev].back().norminf(0,1,IntVect(0));
300 m_is_singular[alev] = (asum <= amax * RT(1.e-12));
301 }
302 }
303 }
304 }
305}
306
307template <typename MF>
308void
309MLALaplacianT<MF>::prepareForSolve ()
310{
311 BL_PROFILE("MLALaplacian::prepareForSolve()");
312 MLCellABecLapT<MF>::prepareForSolve();
313 averageDownCoeffs();
314 updateSingularFlag();
315 m_needs_update = false;
316}
317
318template <typename MF>
319void
320MLALaplacianT<MF>::update ()
321{
322 if (MLCellABecLapT<MF>::needsUpdate()) { MLCellABecLapT<MF>::update(); }
323 averageDownCoeffs();
324 updateSingularFlag();
325 m_needs_update = false;
326}
327
328template <typename MF>
329void
330MLALaplacianT<MF>::Fapply (int amrlev, int mglev, MF& out, const MF& in) const
331{
332 BL_PROFILE("MLALaplacian::Fapply()");
333
334 const int ncomp = this->getNComp();
335
336 const MF& acoef = m_a_coeffs[amrlev][mglev];
337
338 const GpuArray<RT,AMREX_SPACEDIM> dxinv
339 {AMREX_D_DECL(RT(this->m_geom[amrlev][mglev].InvCellSize(0)),
340 RT(this->m_geom[amrlev][mglev].InvCellSize(1)),
341 RT(this->m_geom[amrlev][mglev].InvCellSize(2)))};
342#if (AMREX_SPACEDIM < 3)
343 const RT dx = RT(this->m_geom[amrlev][mglev].CellSize(0));
344 const RT probxlo = RT(this->m_geom[amrlev][mglev].ProbLo(0));
345#endif
346
347#if (AMREX_SPACEDIM == 3)
348 GpuArray<RT,2> dhinv {this->get_d0(dxinv[0], dxinv[1], dxinv[2]),
349 this->get_d1(dxinv[0], dxinv[1], dxinv[2])};
350#endif
351
352 const RT ascalar = m_a_scalar;
353 const RT bscalar = m_b_scalar;
354
355#ifdef AMREX_USE_OMP
356#pragma omp parallel if (Gpu::notInLaunchRegion())
357#endif
358 for (MFIter mfi(out, TilingIfNotGPU()); mfi.isValid(); ++mfi)
359 {
360 const Box& bx = mfi.tilebox();
361 const auto& xfab = in.array(mfi);
362 const auto& yfab = out.array(mfi);
363 const auto& afab = acoef.array(mfi);
364
365#if (AMREX_SPACEDIM != 3)
366 if (this->m_has_metric_term) {
367 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( bx, tbx,
368 {
369 mlalap_adotx_m(tbx, yfab, xfab, afab, dxinv, ascalar, bscalar, dx, probxlo, ncomp);
370 });
371 } else {
372 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( bx, tbx,
373 {
374 mlalap_adotx(tbx, yfab, xfab, afab, dxinv, ascalar, bscalar, ncomp);
375 });
376 }
377#else
378 if (this->hasHiddenDimension()) {
379 Box const& bx2d = this->compactify(bx);
380 const auto& xfab2d = this->compactify(xfab);
381 const auto& yfab2d = this->compactify(yfab);
382 const auto& afab2d = this->compactify(afab);
383 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( bx2d, tbx2d,
384 {
385 TwoD::mlalap_adotx(tbx2d, yfab2d, xfab2d, afab2d, dhinv, ascalar, bscalar, ncomp);
386 });
387 } else {
388 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( bx, tbx,
389 {
390 mlalap_adotx(tbx, yfab, xfab, afab, dxinv, ascalar, bscalar, ncomp);
391 });
392 }
393#endif
394 }
395}
396
397template <typename MF>
398void
399MLALaplacianT<MF>::normalize (int amrlev, int mglev, MF& mf) const
400{
401 BL_PROFILE("MLALaplacian::normalize()");
402
403 const int ncomp = this->getNComp();
404
405 const MF& acoef = m_a_coeffs[amrlev][mglev];
406
407 const GpuArray<RT,AMREX_SPACEDIM> dxinv
408 {AMREX_D_DECL(RT(this->m_geom[amrlev][mglev].InvCellSize(0)),
409 RT(this->m_geom[amrlev][mglev].InvCellSize(1)),
410 RT(this->m_geom[amrlev][mglev].InvCellSize(2)))};
411#if (AMREX_SPACEDIM < 3)
412 const RT dx = RT(this->m_geom[amrlev][mglev].CellSize(0));
413 const RT probxlo = RT(this->m_geom[amrlev][mglev].ProbLo(0));
414#endif
415
416#if (AMREX_SPACEDIM == 3)
417 GpuArray<RT,2> dhinv {this->get_d0(dxinv[0], dxinv[1], dxinv[2]),
418 this->get_d1(dxinv[0], dxinv[1], dxinv[2])};
419#endif
420
421 const RT ascalar = m_a_scalar;
422 const RT bscalar = m_b_scalar;
423
424#ifdef AMREX_USE_OMP
425#pragma omp parallel if (Gpu::notInLaunchRegion())
426#endif
427 for (MFIter mfi(mf, TilingIfNotGPU()); mfi.isValid(); ++mfi)
428 {
429 const Box& bx = mfi.tilebox();
430 const auto& fab = mf.array(mfi);
431 const auto& afab = acoef.array(mfi);
432
433#if (AMREX_SPACEDIM != 3)
434 if (this->m_has_metric_term) {
435 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( bx, tbx,
436 {
437 mlalap_normalize_m(tbx, fab, afab, dxinv, ascalar, bscalar, dx, probxlo, ncomp);
438 });
439 } else {
440 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( bx, tbx,
441 {
442 mlalap_normalize(tbx, fab, afab, dxinv, ascalar, bscalar, ncomp);
443 });
444 }
445#else
446 if (this->hasHiddenDimension()) {
447 Box const& bx2d = this->compactify(bx);
448 const auto& fab2d = this->compactify(fab);
449 const auto& afab2d = this->compactify(afab);
450 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( bx2d, tbx2d,
451 {
452 TwoD::mlalap_normalize(tbx2d, fab2d, afab2d, dhinv, ascalar, bscalar, ncomp);
453 });
454 } else {
455 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( bx, tbx,
456 {
457 mlalap_normalize(tbx, fab, afab, dxinv, ascalar, bscalar, ncomp);
458 });
459 }
460#endif
461 }
462}
463
464template <typename MF>
465void
466MLALaplacianT<MF>::Fsmooth (int amrlev, int mglev, MF& sol, const MF& rhs, int redblack) const
467{
468 BL_PROFILE("MLALaplacian::Fsmooth()");
469
470 const int ncomp = this->getNComp();
471
472 const MF& acoef = m_a_coeffs[amrlev][mglev];
473 const auto& undrrelxr = this->m_undrrelxr[amrlev][mglev];
474 const auto& maskvals = this->m_maskvals [amrlev][mglev];
475
476 OrientationIter oitr;
477
478 const auto& f0 = undrrelxr[oitr()]; ++oitr;
479 const auto& f1 = undrrelxr[oitr()]; ++oitr;
480#if (AMREX_SPACEDIM > 1)
481 const auto& f2 = undrrelxr[oitr()]; ++oitr;
482 const auto& f3 = undrrelxr[oitr()]; ++oitr;
483#if (AMREX_SPACEDIM > 2)
484 const auto& f4 = undrrelxr[oitr()]; ++oitr;
485 const auto& f5 = undrrelxr[oitr()]; ++oitr;
486#endif
487#endif
488
489 const MultiMask& mm0 = maskvals[0];
490 const MultiMask& mm1 = maskvals[1];
491#if (AMREX_SPACEDIM > 1)
492 const MultiMask& mm2 = maskvals[2];
493 const MultiMask& mm3 = maskvals[3];
494#if (AMREX_SPACEDIM > 2)
495 const MultiMask& mm4 = maskvals[4];
496 const MultiMask& mm5 = maskvals[5];
497#endif
498#endif
499
500 const Real* dxinv = this->m_geom[amrlev][mglev].InvCellSize();
501 AMREX_D_TERM(const RT dhx = m_b_scalar*RT(dxinv[0]*dxinv[0]);,
502 const RT dhy = m_b_scalar*RT(dxinv[1]*dxinv[1]);,
503 const RT dhz = m_b_scalar*RT(dxinv[2]*dxinv[2]););
504
505#if (AMREX_SPACEDIM == 3)
506 RT dh0 = this->get_d0(dhx, dhy, dhz);
507 RT dh1 = this->get_d1(dhx, dhy, dhz);
508#endif
509
510#if (AMREX_SPACEDIM < 3)
511 const RT dx = RT(this->m_geom[amrlev][mglev].CellSize(0));
512 const RT probxlo = RT(this->m_geom[amrlev][mglev].ProbLo(0));
513#endif
514
515 const RT alpha = m_a_scalar;
516
517 MFItInfo mfi_info;
518 if (Gpu::notInLaunchRegion()) { mfi_info.EnableTiling().SetDynamic(true); }
519
520#ifdef AMREX_USE_OMP
521#pragma omp parallel if (Gpu::notInLaunchRegion())
522#endif
523 for (MFIter mfi(sol,mfi_info); mfi.isValid(); ++mfi)
524 {
525 const auto& m0 = mm0.array(mfi);
526 const auto& m1 = mm1.array(mfi);
527#if (AMREX_SPACEDIM > 1)
528 const auto& m2 = mm2.array(mfi);
529 const auto& m3 = mm3.array(mfi);
530#if (AMREX_SPACEDIM > 2)
531 const auto& m4 = mm4.array(mfi);
532 const auto& m5 = mm5.array(mfi);
533#endif
534#endif
535
536 const Box& tbx = mfi.tilebox();
537 const Box& vbx = mfi.validbox();
538 const auto& solnfab = sol.array(mfi);
539 const auto& rhsfab = rhs.array(mfi);
540 const auto& afab = acoef.array(mfi);
541
542 const auto& f0fab = f0.array(mfi);
543 const auto& f1fab = f1.array(mfi);
544#if (AMREX_SPACEDIM > 1)
545 const auto& f2fab = f2.array(mfi);
546 const auto& f3fab = f3.array(mfi);
547#if (AMREX_SPACEDIM > 2)
548 const auto& f4fab = f4.array(mfi);
549 const auto& f5fab = f5.array(mfi);
550#endif
551#endif
552
553#if (AMREX_SPACEDIM == 1)
554 if (this->m_has_metric_term) {
555 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( tbx, thread_box,
556 {
557 mlalap_gsrb_m(thread_box, solnfab, rhsfab, alpha, dhx,
558 afab,
559 f0fab, m0,
560 f1fab, m1,
561 vbx, redblack,
562 dx, probxlo, ncomp);
563 });
564 } else {
565 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( tbx, thread_box,
566 {
567 mlalap_gsrb(thread_box, solnfab, rhsfab, alpha, dhx,
568 afab,
569 f0fab, m0,
570 f1fab, m1,
571 vbx, redblack, ncomp);
572 });
573 }
574
575#endif
576
577#if (AMREX_SPACEDIM == 2)
578 if (this->m_has_metric_term) {
579 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( tbx, thread_box,
580 {
581 mlalap_gsrb_m(thread_box, solnfab, rhsfab, alpha, dhx, dhy,
582 afab,
583 f0fab, m0,
584 f1fab, m1,
585 f2fab, m2,
586 f3fab, m3,
587 vbx, redblack,
588 dx, probxlo, ncomp);
589 });
590 } else {
591 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( tbx, thread_box,
592 {
593 mlalap_gsrb(thread_box, solnfab, rhsfab, alpha, dhx, dhy,
594 afab,
595 f0fab, m0,
596 f1fab, m1,
597 f2fab, m2,
598 f3fab, m3,
599 vbx, redblack, ncomp);
600 });
601 }
602#endif
603
604#if (AMREX_SPACEDIM == 3)
605 if (this->hasHiddenDimension()) {
606 Box const& tbx_2d = this->compactify(tbx);
607 Box const& vbx_2d = this->compactify(vbx);
608 const auto& solnfab_2d = this->compactify(solnfab);
609 const auto& rhsfab_2d = this->compactify(rhsfab);
610 const auto& afab_2d = this->compactify(afab);
611 const auto& f0fab_2d = this->compactify(this->get_d0(f0fab,f1fab,f2fab));
612 const auto& f1fab_2d = this->compactify(this->get_d1(f0fab,f1fab,f2fab));
613 const auto& f2fab_2d = this->compactify(this->get_d0(f3fab,f4fab,f5fab));
614 const auto& f3fab_2d = this->compactify(this->get_d1(f3fab,f4fab,f5fab));
615 const auto& m0_2d = this->compactify(this->get_d0(m0,m1,m2));
616 const auto& m1_2d = this->compactify(this->get_d1(m0,m1,m2));
617 const auto& m2_2d = this->compactify(this->get_d0(m3,m4,m5));
618 const auto& m3_2d = this->compactify(this->get_d1(m3,m4,m5));
619 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( tbx_2d, thread_box,
620 {
621 TwoD::mlalap_gsrb(thread_box, solnfab_2d, rhsfab_2d, alpha, dh0, dh1,
622 afab_2d,
623 f0fab_2d, m0_2d,
624 f1fab_2d, m1_2d,
625 f2fab_2d, m2_2d,
626 f3fab_2d, m3_2d,
627 vbx_2d, redblack, ncomp);
628 });
629 } else {
630 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( tbx, thread_box,
631 {
632 mlalap_gsrb(thread_box, solnfab, rhsfab, alpha, dhx, dhy, dhz,
633 afab,
634 f0fab, m0,
635 f1fab, m1,
636 f2fab, m2,
637 f3fab, m3,
638 f4fab, m4,
639 f5fab, m5,
640 vbx, redblack, ncomp);
641 });
642 }
643#endif
644 }
645}
646
647template <typename MF>
648void
649MLALaplacianT<MF>::FFlux (int amrlev, const MFIter& mfi,
650 const Array<FAB*,AMREX_SPACEDIM>& flux,
651 const FAB& sol, Location, int face_only) const
652{
653 BL_PROFILE("MLALaplacian::FFlux()");
654
655 const int ncomp = this->getNComp();
656 const int mglev = 0;
657 const Box& box = mfi.tilebox();
658 const Real* dxinv = this->m_geom[amrlev][mglev].InvCellSize();
659
660 AMREX_D_TERM(const auto& fxarr = flux[0]->array();,
661 const auto& fyarr = flux[1]->array();,
662 const auto& fzarr = flux[2]->array(););
663 const auto& solarr = sol.array();
664
665#if (AMREX_SPACEDIM != 3)
666 const RT dx = RT(this->m_geom[amrlev][mglev].CellSize(0));
667 const RT probxlo = RT(this->m_geom[amrlev][mglev].ProbLo(0));
668#endif
669
670#if (AMREX_SPACEDIM == 3)
671 if (face_only) {
672 if (this->hiddenDirection() != 0) {
673 RT fac = m_b_scalar * RT(dxinv[0]);
674 Box blo = amrex::bdryLo(box, 0);
675 int blen = box.length(0);
676 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( blo, tbox,
677 {
678 mlalap_flux_xface(tbox, fxarr, solarr, fac, blen, ncomp);
679 });
680 } else {
681 flux[0]->template setVal<RunOn::Device>(RT(0.0));
682 }
683 if (this->hiddenDirection() != 1) {
684 RT fac = m_b_scalar * RT(dxinv[1]);
685 Box blo = amrex::bdryLo(box, 1);
686 int blen = box.length(1);
687 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( blo, tbox,
688 {
689 mlalap_flux_yface(tbox, fyarr, solarr, fac, blen, ncomp);
690 });
691 } else {
692 flux[1]->template setVal<RunOn::Device>(RT(0.0));
693 }
694 if (this->hiddenDirection() != 2) {
695 RT fac = m_b_scalar * RT(dxinv[2]);
696 Box blo = amrex::bdryLo(box, 2);
697 int blen = box.length(2);
698 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( blo, tbox,
699 {
700 mlalap_flux_zface(tbox, fzarr, solarr, fac, blen, ncomp);
701 });
702 } else {
703 flux[2]->template setVal<RunOn::Device>(RT(0.0));
704 }
705 } else {
706 if (this->hiddenDirection() != 0) {
707 RT fac = m_b_scalar * RT(dxinv[0]);
708 Box bflux = amrex::surroundingNodes(box, 0);
709 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( bflux, tbox,
710 {
711 mlalap_flux_x(tbox, fxarr, solarr, fac, ncomp);
712 });
713 } else {
714 flux[0]->template setVal<RunOn::Device>(RT(0.0));
715 }
716 if (this->hiddenDirection() != 1) {
717 RT fac = m_b_scalar * RT(dxinv[1]);
718 Box bflux = amrex::surroundingNodes(box, 1);
719 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( bflux, tbox,
720 {
721 mlalap_flux_y(tbox, fyarr, solarr, fac, ncomp);
722 });
723 } else {
724 flux[1]->template setVal<RunOn::Device>(RT(0.0));
725 }
726 if (this->hiddenDirection() != 2) {
727 RT fac = m_b_scalar * RT(dxinv[2]);
728 Box bflux = amrex::surroundingNodes(box, 2);
729 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( bflux, tbox,
730 {
731 mlalap_flux_z(tbox, fzarr, solarr, fac, ncomp);
732 });
733 } else {
734 flux[2]->template setVal<RunOn::Device>(RT(0.0));
735 }
736 }
737#elif (AMREX_SPACEDIM == 2)
738 if (face_only) {
739 if (this->hiddenDirection() != 0) {
740 RT fac = m_b_scalar * RT(dxinv[0]);
741 Box blo = amrex::bdryLo(box, 0);
742 int blen = box.length(0);
743 if (this->m_has_metric_term) {
744 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( blo, tbox,
745 {
746 mlalap_flux_xface_m(tbox, fxarr, solarr, fac, blen, dx, probxlo, ncomp);
747 });
748 } else {
749 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( blo, tbox,
750 {
751 mlalap_flux_xface(tbox, fxarr, solarr, fac, blen, ncomp);
752 });
753 }
754 } else {
755 flux[0]->template setVal<RunOn::Device>(RT(0.0));
756 }
757 if (this->hiddenDirection() != 1) {
758 RT fac = m_b_scalar * RT(dxinv[1]);
759 Box blo = amrex::bdryLo(box, 1);
760 int blen = box.length(1);
761 if (this->m_has_metric_term) {
762 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( blo, tbox,
763 {
764 mlalap_flux_yface_m(tbox, fyarr, solarr, fac, blen, dx, probxlo, ncomp);
765 });
766 } else {
767 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( blo, tbox,
768 {
769 mlalap_flux_yface(tbox, fyarr, solarr, fac, blen, ncomp);
770 });
771 }
772 } else {
773 flux[1]->template setVal<RunOn::Device>(RT(0.0));
774 }
775 } else {
776 if (this->hiddenDirection() != 0) {
777 RT fac = m_b_scalar * RT(dxinv[0]);
778 Box bflux = amrex::surroundingNodes(box, 0);
779 if (this->m_has_metric_term) {
780 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( bflux, tbox,
781 {
782 mlalap_flux_x_m(tbox, fxarr, solarr, fac, dx, probxlo, ncomp);
783 });
784 } else {
785 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( bflux, tbox,
786 {
787 mlalap_flux_x(tbox, fxarr, solarr, fac, ncomp);
788 });
789 }
790 } else {
791 flux[0]->template setVal<RunOn::Device>(RT(0.0));
792 }
793 if (this->hiddenDirection() != 1) {
794 RT fac = m_b_scalar * RT(dxinv[1]);
795 Box bflux = amrex::surroundingNodes(box, 1);
796 if (this->m_has_metric_term) {
797 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( bflux, tbox,
798 {
799 mlalap_flux_y_m(tbox, fyarr, solarr, fac, dx, probxlo, ncomp);
800 });
801 } else {
802 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( bflux, tbox,
803 {
804 mlalap_flux_y(tbox, fyarr, solarr, fac, ncomp);
805 });
806 }
807 } else {
808 flux[1]->template setVal<RunOn::Device>(RT(0.0));
809 }
810 }
811#else
812 if (face_only) {
813 RT fac = m_b_scalar * RT(dxinv[0]);
814 Box blo = amrex::bdryLo(box, 0);
815 int blen = box.length(0);
816 if (this->m_has_metric_term) {
817 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( blo, tbox,
818 {
819 mlalap_flux_xface_m(tbox, fxarr, solarr, fac, blen, dx, probxlo, ncomp);
820 });
821 } else {
822 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( blo, tbox,
823 {
824 mlalap_flux_xface(tbox, fxarr, solarr, fac, blen, ncomp);
825 });
826 }
827 } else {
828 RT fac = m_b_scalar * RT(dxinv[0]);
829 Box bflux = amrex::surroundingNodes(box, 0);
830 if (this->m_has_metric_term) {
831 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( bflux, tbox,
832 {
833 mlalap_flux_x_m(tbox, fxarr, solarr, fac, dx, probxlo, ncomp);
834 });
835 } else {
836 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE ( bflux, tbox,
837 {
838 mlalap_flux_x(tbox, fxarr, solarr, fac, ncomp);
839 });
840 }
841 }
842#endif
843}
844
845extern template class MLALaplacianT<MultiFab>;
846
847using MLALaplacian = MLALaplacianT<MultiFab>;
848
849}
850
851#endif
BL_PROFILE
#define BL_PROFILE(a)
Definition AMReX_BLProfiler.H:551
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
Definition AMReX_FabFactory.H:50
amrex::MFIter
Iterator for looping ever tiles and boxes of amrex::FabArray based containers.
Definition AMReX_MFIter.H:88
amrex::MFIter::tilebox
Box tilebox() const noexcept
Return the tile Box at the current index.
Definition AMReX_MFIter.cpp:389
amrex::MFIter::isValid
bool isValid() const noexcept
Is the iterator valid i.e. is it associated with a FAB?
Definition AMReX_MFIter.H:172
amrex::MLALaplacianT
Multi-component ALaplacian (a scalar plus optional spatial a coeffs).
Definition AMReX_MLALaplacian.H:22
amrex::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:263
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:649
amrex::MLALaplacianT::update
void update() override
Update for reuse.
Definition AMReX_MLALaplacian.H:320
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:309
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:330
amrex::MLALaplacianT::averageDownCoeffs
void averageDownCoeffs()
Average a coefficients down across all AMR and MG levels.
Definition AMReX_MLALaplacian.H:225
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:399
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:466
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:242
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:216
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:305
amrex::MLCellABecLapT::update
void update() override
Average coefficients/metrics when marked dirty.
Definition AMReX_MLCellABecLap.H:298
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:79
amrex::surroundingNodes
__host__ __device__ BoxND< dim > surroundingNodes(const BoxND< dim > &b, int dir) noexcept
Return a BoxND with NODE based coordinates in direction dir that encloses BoxND b.
Definition AMReX_Box.H:1582
amrex::bdryLo
__host__ __device__ BoxND< dim > bdryLo(const BoxND< dim > &b, int dir, int len=1) noexcept
Return the BoxND of length len on the low boundary of b along coordinate direction dir.
Definition AMReX_Box.H:1715
amrex::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:241
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