Block-Structured AMR Software Framework
AMReX_MLEBNodeFDLaplacian.H
Go to the documentation of this file.
1 #ifndef AMREX_MLEBNODEFDLAPLACIAN_H_
2 #define AMREX_MLEBNODEFDLAPLACIAN_H_
3 #include <AMReX_Config.H>
4 
5 #include <AMReX_Array.H>
6 #ifdef AMREX_USE_EB
7 #include <AMReX_EBFabFactory.H>
8 #endif
9 #include <AMReX_MLNodeLinOp.H>
10 
11 #include <limits>
12 
13 namespace amrex {
14 
15 // Although the class has EB in the name, it works for non-EB build too.
16 //
17 // del dot (sigma grad phi) = rhs, for non-RZ
18 // where phi and rhs are nodal multifab, and sigma is a tensor constant
19 // with only diagonal components. The EB is assumed to be Dirichlet.
20 //
21 // del dot (sigma grad phi) - alpha/r^2 phi = rhs, for RZ where alpha is a
22 // scalar constant that is zero by default.
23 //
24 // New feature: for non-RZ, sigma can also be a single-component
25 // cell-centered multifab.
26 
28  : public MLNodeLinOp
29 {
30 public:
31 
32  MLEBNodeFDLaplacian () = default;
33 
34 #ifdef AMREX_USE_EB
36  const Vector<BoxArray>& a_grids,
37  const Vector<DistributionMapping>& a_dmap,
38  const LPInfo& a_info,
39  const Vector<EBFArrayBoxFactory const*>& a_factory);
40 #endif
41 
43  const Vector<BoxArray>& a_grids,
44  const Vector<DistributionMapping>& a_dmap,
45  const LPInfo& a_info);
46 
47  ~MLEBNodeFDLaplacian () override = default;
48 
53 
54  void setSigma (Array<Real,AMREX_SPACEDIM> const& a_sigma) noexcept;
55 
56  void setSigma (int amrlev, MultiFab const& a_sigma);
57 
58  void setRZ (bool flag);
59 
60  void setAlpha (Real a_alpha);
61 
62 #ifdef AMREX_USE_EB
63 
64  // Phi on EB
65  void setEBDirichlet (Real a_phi_eb);
66  //
67  template <typename F>
68  std::enable_if_t<IsCallableR<Real,F,AMREX_D_DECL(Real,Real,Real)>::value>
69  setEBDirichlet (F const& f);
70 
71  void define (const Vector<Geometry>& a_geom,
72  const Vector<BoxArray>& a_grids,
73  const Vector<DistributionMapping>& a_dmap,
74  const LPInfo& a_info,
75  const Vector<EBFArrayBoxFactory const*>& a_factory);
76 
77  [[nodiscard]] std::unique_ptr<FabFactory<FArrayBox> > makeFactory (int amrlev, int mglev) const final;
78 
79  void scaleRHS (int amrlev, MultiFab& rhs) const final;
80 
81 #endif
82 
83  void define (const Vector<Geometry>& a_geom,
84  const Vector<BoxArray>& a_grids,
85  const Vector<DistributionMapping>& a_dmap,
86  const LPInfo& a_info);
87 
88  [[nodiscard]] std::string name () const override { return std::string("MLEBNodeFDLaplacian"); }
89 
90  void restriction (int amrlev, int cmglev, MultiFab& crse, MultiFab& fine) const final;
91  void interpolation (int amrlev, int fmglev, MultiFab& fine, const MultiFab& crse) const final;
92 
93  void prepareForSolve () final;
94  void Fapply (int amrlev, int mglev, MultiFab& out, const MultiFab& in) const final;
95  void Fsmooth (int amrlev, int mglev, MultiFab& sol, const MultiFab& rhs) const final;
96  void normalize (int amrlev, int mglev, MultiFab& mf) const final;
97 
98  void fixUpResidualMask (int amrlev, iMultiFab& resmsk) final;
99 
100  [[nodiscard]] bool isSingular (int) const final { return false; }
101  [[nodiscard]] bool isBottomSingular () const final { return false; }
102 
103  void compGrad (int amrlev, const Array<MultiFab*,AMREX_SPACEDIM>& grad,
104  MultiFab& sol, Location /*loc*/) const override;
105 
106 #if defined(AMREX_USE_HYPRE) && (AMREX_SPACEDIM > 1)
107  void fillIJMatrix (MFIter const& mfi,
109  Array4<int const> const& lid,
110  HypreNodeLap::Int* ncols,
111  HypreNodeLap::Int* cols,
112  Real* mat) const override;
113 
114  void fillRHS (MFIter const& mfi,
115  Array4<int const> const& lid,
116  Real* rhs,
117  Array4<Real const> const& bfab) const override;
118 #endif
119 
120  void postSolve (Vector<MultiFab>& sol) const override;
121 
122 private:
125  bool m_has_sigma_mf = false;
126  Real m_s_phi_eb = std::numeric_limits<Real>::lowest();
128  int m_rz = false;
129  Real m_rz_alpha = 0._rt;
130 };
131 
132 #ifdef AMREX_USE_EB
133 
134 template <typename F>
135 std::enable_if_t<IsCallableR<Real,F,AMREX_D_DECL(Real,Real,Real)>::value>
136 MLEBNodeFDLaplacian::setEBDirichlet (F const& f)
137 {
138  m_phi_eb.resize(m_num_amr_levels);
139  for (int amrlev = 0; amrlev < m_num_amr_levels; ++amrlev) {
140  auto const* factory = dynamic_cast<EBFArrayBoxFactory const*>(m_factory[amrlev][0].get());
141  if (factory) {
142  Geometry const& geom = m_geom[amrlev][0];
143  auto const problo = geom.ProbLoArray();
144  auto const cellsize = geom.CellSizeArray();
145  m_phi_eb[amrlev].define(amrex::convert(m_grids[amrlev][0],IntVect(1)),
146  m_dmap[amrlev][0], 1, 1);
147  m_phi_eb[amrlev].setVal(0.0);
148  auto const& flags = factory->getMultiEBCellFlagFab();
149  auto const& levset = factory->getLevelSet();
150 #ifdef AMREX_USE_OMP
151 #pragma omp parallel if (Gpu::notInLaunchRegion())
152 #endif
153  for (MFIter mfi(m_phi_eb[amrlev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
154  {
155  const Box& ndbx = mfi.growntilebox();
156  const auto& flag = flags[mfi];
157  if (flag.getType() != FabType::regular) {
158  Array4<Real const> const lstarr = levset.const_array(mfi);
159  Array4<Real> const& phi = m_phi_eb[amrlev].array(mfi);
160  AMREX_HOST_DEVICE_FOR_3D(ndbx, i, j, k,
161  {
162  if (lstarr(i,j,k) >= Real(0.0)) {
163  phi(i,j,k) = f(AMREX_D_DECL(problo[0]+i*cellsize[0],
164  problo[1]+j*cellsize[1],
165  problo[2]+k*cellsize[2]));
166  }
167  });
168  }
169  }
170  }
171  }
172 }
173 
174 #endif
175 
176 }
177 
178 #endif
#define AMREX_HOST_DEVICE_FOR_3D(...)
Definition: AMReX_GpuLaunch.nolint.H:50
Array4< Real > fine
Definition: AMReX_InterpFaceRegister.cpp:90
Array4< Real const > crse
Definition: AMReX_InterpFaceRegister.cpp:92
#define AMREX_D_DECL(a, b, c)
Definition: AMReX_SPACE.H:104
GpuArray< Real, AMREX_SPACEDIM > CellSizeArray() const noexcept
Definition: AMReX_CoordSys.H:76
Definition: AMReX_EBFabFactory.H:22
Rectangular problem domain geometry.
Definition: AMReX_Geometry.H:73
GpuArray< Real, AMREX_SPACEDIM > ProbLoArray() const noexcept
Definition: AMReX_Geometry.H:186
HYPRE_Int Int
Definition: AMReX_HypreNodeLap.H:36
Definition: AMReX_MFIter.H:57
bool isValid() const noexcept
Is the iterator valid i.e. is it associated with a FAB?
Definition: AMReX_MFIter.H:141
Definition: AMReX_MLEBNodeFDLaplacian.H:29
void restriction(int amrlev, int cmglev, MultiFab &crse, MultiFab &fine) const final
Definition: AMReX_MLEBNodeFDLaplacian.cpp:167
std::string name() const override
Definition: AMReX_MLEBNodeFDLaplacian.H:88
bool isBottomSingular() const final
Is the bottom of MG singular?
Definition: AMReX_MLEBNodeFDLaplacian.H:101
void postSolve(Vector< MultiFab > &sol) const override
Definition: AMReX_MLEBNodeFDLaplacian.cpp:754
int m_rz
Definition: AMReX_MLEBNodeFDLaplacian.H:128
void Fsmooth(int amrlev, int mglev, MultiFab &sol, const MultiFab &rhs) const final
Definition: AMReX_MLEBNodeFDLaplacian.cpp:521
void interpolation(int amrlev, int fmglev, MultiFab &fine, const MultiFab &crse) const final
Definition: AMReX_MLEBNodeFDLaplacian.cpp:214
void setAlpha(Real a_alpha)
Definition: AMReX_MLEBNodeFDLaplacian.cpp:67
MLEBNodeFDLaplacian(const MLEBNodeFDLaplacian &)=delete
Real m_rz_alpha
Definition: AMReX_MLEBNodeFDLaplacian.H:129
void setRZ(bool flag)
Definition: AMReX_MLEBNodeFDLaplacian.cpp:57
GpuArray< Real, AMREX_SPACEDIM > m_sigma
Definition: AMReX_MLEBNodeFDLaplacian.H:123
void Fapply(int amrlev, int mglev, MultiFab &out, const MultiFab &in) const final
Definition: AMReX_MLEBNodeFDLaplacian.cpp:398
void define(const Vector< Geometry > &a_geom, const Vector< BoxArray > &a_grids, const Vector< DistributionMapping > &a_dmap, const LPInfo &a_info)
Definition: AMReX_MLEBNodeFDLaplacian.cpp:123
MLEBNodeFDLaplacian & operator=(const MLEBNodeFDLaplacian &)=delete
Vector< MultiFab > m_phi_eb
Definition: AMReX_MLEBNodeFDLaplacian.H:127
void setSigma(Array< Real, AMREX_SPACEDIM > const &a_sigma) noexcept
Definition: AMReX_MLEBNodeFDLaplacian.cpp:36
void normalize(int amrlev, int mglev, MultiFab &mf) const final
Definition: AMReX_MLEBNodeFDLaplacian.cpp:626
MLEBNodeFDLaplacian(MLEBNodeFDLaplacian &&)=delete
Real m_s_phi_eb
Definition: AMReX_MLEBNodeFDLaplacian.H:126
void prepareForSolve() final
Definition: AMReX_MLEBNodeFDLaplacian.cpp:258
void compGrad(int amrlev, const Array< MultiFab *, AMREX_SPACEDIM > &grad, MultiFab &sol, Location) const override
Definition: AMReX_MLEBNodeFDLaplacian.cpp:638
~MLEBNodeFDLaplacian() override=default
bool m_has_sigma_mf
Definition: AMReX_MLEBNodeFDLaplacian.H:125
Vector< Vector< std::unique_ptr< MultiFab > > > m_sigma_mf
Definition: AMReX_MLEBNodeFDLaplacian.H:124
void fixUpResidualMask(int amrlev, iMultiFab &resmsk) final
Definition: AMReX_MLEBNodeFDLaplacian.cpp:632
bool isSingular(int) const final
Is it singular on given AMR level?
Definition: AMReX_MLEBNodeFDLaplacian.H:100
Vector< Vector< std::unique_ptr< FabFactory< FAB > > > > m_factory
Definition: AMReX_MLLinOp.H:601
Vector< Vector< BoxArray > > m_grids
Definition: AMReX_MLLinOp.H:599
virtual std::unique_ptr< FabFactory< FAB > > makeFactory(int, int) const
Definition: AMReX_MLLinOp.H:689
Vector< Vector< DistributionMapping > > m_dmap
Definition: AMReX_MLLinOp.H:600
Vector< Vector< Geometry > > m_geom
first Vector is for amr level and second is mg level
Definition: AMReX_MLLinOp.H:598
LinOpEnumType::Location Location
Definition: AMReX_MLLinOp.H:114
virtual void scaleRHS(int, MF &) const
scale RHS to fix solvability
Definition: AMReX_MLLinOp.H:458
int m_num_amr_levels
Definition: AMReX_MLLinOp.H:583
Definition: AMReX_MLNodeLinOp.H:16
A collection (stored as an array) of FArrayBox objects.
Definition: AMReX_MultiFab.H:38
This class is a thin wrapper around std::vector. Unlike vector, Vector::operator[] provides bound che...
Definition: AMReX_Vector.H:27
Definition: AMReX_iMultiFab.H:32
static int f(amrex::Real t, N_Vector y_data, N_Vector y_rhs, void *user_data)
Definition: AMReX_SundialsIntegrator.H:44
Definition: AMReX_Amr.cpp:49
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE BoxND< dim > convert(const BoxND< dim > &b, const IntVectND< dim > &typ) noexcept
Returns a BoxND with different type.
Definition: AMReX_Box.H:1435
IntVectND< AMREX_SPACEDIM > IntVect
Definition: AMReX_BaseFwd.H:30
bool TilingIfNotGPU() noexcept
Definition: AMReX_MFIter.H:12
std::array< T, N > Array
Definition: AMReX_Array.H:24
Definition: AMReX_Array4.H:61
Test if a given type T is callable with arguments of type Args...
Definition: AMReX_TypeTraits.H:207
Definition: AMReX_MLLinOp.H:35