ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
ERF_PlaneAverage.H
Go to the documentation of this file.
1 #ifndef ERF_PlaneAverage_H
2 #define ERF_PlaneAverage_H
3 
4 #include "AMReX_Gpu.H"
5 #include "AMReX_iMultiFab.H"
6 #include "AMReX_MultiFab.H"
7 #include "AMReX_GpuContainers.H"
9 
10 /**
11  * Basic averaging and interpolation operations
12  */
13 
14 class PlaneAverage {
15 public:
16  AMREX_FORCE_INLINE
17  explicit PlaneAverage (const amrex::MultiFab* field_in,
18  amrex::Geometry geom_in,
19  int axis_in,
20  amrex::IntVect n_ghost_to_inc=amrex::IntVect{0});
21  PlaneAverage () = delete;
22  ~PlaneAverage () = default;
23 
24  AMREX_FORCE_INLINE
25  void operator()();
26 
27  /** evaluate line average at specific location for any average component */
28  [[nodiscard]] AMREX_FORCE_INLINE
29  amrex::Real line_average_interpolated (amrex::Real x, int comp) const;
30 
31  /** change precision of text file output */
32  void set_precision (int p) { m_precision = p; }
33 
34  [[nodiscard]] amrex::Real dx () const { return m_dx; }
35  [[nodiscard]] amrex::Real xlo () const { return m_xlo; }
36 
37  [[nodiscard]] int axis () const { return m_axis; }
38  [[nodiscard]] int level () const { return m_level; }
39  [[nodiscard]] int ncomp () const { return m_ncomp; }
40  [[nodiscard]] int ncell_plane () const { return m_ncell_plane; }
41  [[nodiscard]] int ncell_line () const { return m_ncell_line; }
42 
43  [[nodiscard]] const amrex::Vector<amrex::Real>& line_average () const
44  {
45  return m_line_average;
46  }
47 
48  AMREX_FORCE_INLINE
49  void line_average (int comp, amrex::Gpu::HostVector<amrex::Real>& l_vec);
50 
51  [[nodiscard]] const amrex::Vector<amrex::Real>& line_centroids () const
52  {
53  return m_line_xcentroid;
54  }
55 
56  [[nodiscard]] const amrex::MultiFab& field () const { return *m_field; }
57 
58 protected:
59  int m_ncomp; /** number of average components */
60 
61  /** line storage for the average velocity and tracer variables */
62  amrex::Vector<amrex::Real> m_line_average;
63 
64  amrex::Vector<amrex::Real> m_line_xcentroid; /** line storage for centroids of each cell along a line*/
65 
66  amrex::Real m_dx; /** mesh spacing in axis direction*/
67  amrex::Real m_xlo; /** bottom of domain in axis direction */
68 
69  int m_ncell_plane; /** number of cells in plane */
70  int m_ncell_line; /** number of cells along line */
71 
72  int m_precision = 4; /** precision for line plot text file */
73  const int m_level = 0; /** level for plane averaging for now fixed at level=0 */
74 
75  const amrex::MultiFab* m_field;
76  amrex::Geometry m_geom;
77  const int m_axis;
78  amrex::IntVect m_ghost_to_inc = amrex::IntVect{0};
79  amrex::IntVect m_ixtype;
80 
81 public:
82  /** fill line storage with averages */
83  template <typename IndexSelector>
84  AMREX_FORCE_INLINE
85  void compute_averages (const IndexSelector& idxOp, const amrex::MultiFab& mfab);
86 };
87 
88 
89 PlaneAverage::PlaneAverage (const amrex::MultiFab* field_in,
90  amrex::Geometry geom_in,
91  int axis_in,
92  amrex::IntVect n_ghost_to_inc)
93  : m_field(field_in), m_geom(geom_in), m_axis(axis_in), m_ghost_to_inc(n_ghost_to_inc)
94 {
95  AMREX_ALWAYS_ASSERT(m_axis >= 0 && m_axis < AMREX_SPACEDIM);
96 
97  m_xlo = m_geom.ProbLo(m_axis);
98  m_dx = m_geom.CellSize(m_axis);
99  m_ncomp = m_field->nComp();
100  m_ixtype = m_field->boxArray().ixType().toIntVect();
101 
102  amrex::Box domain = m_geom.Domain();
103  m_ghost_to_inc = std::min(field_in->nGrowVect(),m_ghost_to_inc);
104  domain.grow(axis_in, m_ghost_to_inc[axis_in]);
105 
106  amrex::IntVect dom_lo(domain.loVect());
107  amrex::IntVect dom_hi(domain.hiVect());
108 
110 
111  m_ncell_plane = 1;
112  auto period = m_geom.periodicity();
113  for (int i = 0; i < AMREX_SPACEDIM; ++i) {
114  int p_fac = (!period.isPeriodic(i)) ? 1 : 0;
115  if (i != m_axis) m_ncell_plane *= (dom_hi[i] - dom_lo[i] + 1 + p_fac*m_ixtype[i]);
116  }
117 
118  m_line_average.resize(static_cast<size_t>(m_ncell_line) * m_ncomp, 0.0);
120 
121  for (int i = 0; i < m_ncell_line; ++i) {
122  m_line_xcentroid[i] = m_xlo + (i + 0.5) * m_dx;
123  }
124 }
125 
126 amrex::Real
127 PlaneAverage::line_average_interpolated (amrex::Real x, int comp) const
128 {
129  AMREX_ALWAYS_ASSERT(comp >= 0 && comp < m_ncomp);
130 
131  amrex::Real c = 0.0;
132  int ind = 0;
133 
134  if (x > m_xlo + 0.5 * m_dx) {
135  ind = static_cast<int>(floor((x - m_xlo) / m_dx - 0.5));
136  const amrex::Real x1 = m_xlo + (ind + 0.5) * m_dx;
137  c = (x - x1) / m_dx;
138  }
139 
140  if (ind + 1 >= m_ncell_line) {
141  ind = m_ncell_line - 2;
142  c = 1.0;
143  }
144 
145  AMREX_ALWAYS_ASSERT(ind >= 0 && ind + 1 < m_ncell_line);
146 
147  return m_line_average[m_ncomp * ind + comp] * (1.0 - c) +
148  m_line_average[m_ncomp * (ind + 1) + comp] * c;
149 }
150 
151 void
152 PlaneAverage::line_average (int comp, amrex::Gpu::HostVector<amrex::Real>& l_vec)
153 {
154  AMREX_ALWAYS_ASSERT(comp >= 0 && comp < m_ncomp);
155 
156  for (int i = 0; i < m_ncell_line; i++)
157  l_vec[i] = m_line_average[m_ncomp * i + comp];
158 }
159 
160 void
162 {
163  std::fill(m_line_average.begin(), m_line_average.end(), 0.0);
164  switch (m_axis) {
165  case 0:
167  break;
168  case 1:
170  break;
171  case 2:
173  break;
174  default:
175  amrex::Abort("axis must be equal to 0, 1, or 2");
176  break;
177  }
178 }
179 
180 template <typename IndexSelector>
181 void
182 PlaneAverage::compute_averages (const IndexSelector& idxOp, const amrex::MultiFab& mfab)
183 {
184  const amrex::Real denom = 1.0 / (amrex::Real)m_ncell_plane;
185  amrex::AsyncArray<amrex::Real> lavg(m_line_average.data(), m_line_average.size());
186  amrex::Real* line_avg = lavg.data();
187  const int ncomp = m_ncomp;
188 
189  amrex::Box domain = amrex::convert(m_geom.Domain(),m_ixtype);
190 
191  amrex::IntVect ng = amrex::IntVect(0);
193  ng[m_axis] = offset;
194 
195  std::unique_ptr<amrex::iMultiFab> mask = OwnerMask(*m_field, m_geom.periodicity(), ng);
196 
197 #ifdef _OPENMP
198 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
199 #endif
200  for (amrex::MFIter mfi(mfab, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) {
201  amrex::Box tbx = mfi.tilebox();
202  if (tbx.smallEnd(m_axis) == domain.smallEnd(m_axis)) tbx.growLo(m_axis,offset);
203  if (tbx.bigEnd (m_axis) == domain.bigEnd (m_axis)) tbx.growHi(m_axis,offset);
204  amrex::Box pbx = PerpendicularBox<IndexSelector>(tbx, amrex::IntVect(0));
205 
206  const amrex::Array4<const amrex::Real>& fab_arr = mfab.const_array(mfi);
207  const amrex::Array4<const int >& mask_arr = mask->const_array(mfi);
208 
209  amrex::ParallelFor(amrex::Gpu::KernelInfo().setReduction(true), pbx, [=]
210  AMREX_GPU_DEVICE( int p_i, int p_j, int p_k,
211  amrex::Gpu::Handler const& handler) noexcept
212  {
213  // Loop over the direction perpendicular to the plane.
214  // This reduces the atomic pressure on the destination arrays.
215  amrex::Box lbx = ParallelBox<IndexSelector>(tbx, amrex::IntVect{p_i, p_j, p_k});
216 
217  for (int k = lbx.smallEnd(2); k <= lbx.bigEnd(2); ++k) {
218  for (int j = lbx.smallEnd(1); j <= lbx.bigEnd(1); ++j) {
219  for (int i = lbx.smallEnd(0); i <= lbx.bigEnd(0); ++i) {
220  int ind = idxOp.getIndx(i, j, k) + offset;
221  for (int n = 0; n < ncomp; ++n) {
222  // NOTE: This factor is to avoid an if statement that will break
223  // the devicereducesum since all threads won't participate.
224  // This more performant than Gpu::Atomic::Add.
225  amrex::Real fac = (mask_arr(i,j,k)) ? 1.0 : 0.0;
226  amrex::Gpu::deviceReduceSum(&line_avg[ncomp * ind + n],
227  fab_arr(i, j, k, n) * denom * fac, handler);
228  }
229  }
230  }
231  }
232  });
233  }
234 
235  lavg.copyToHost(m_line_average.data(), m_line_average.size());
236  amrex::ParallelDescriptor::ReduceRealSum(m_line_average.data(), m_line_average.size());
237 }
238 #endif /* ERF_PlaneAverage.H */
const auto & dom_hi
Definition: ERF_DiffSetup.H:10
const auto & dom_lo
Definition: ERF_DiffSetup.H:9
DirectionSelector< 0 > XDir
Definition: ERF_DirectionSelector.H:36
DirectionSelector< 1 > YDir
Definition: ERF_DirectionSelector.H:37
DirectionSelector< 2 > ZDir
Definition: ERF_DirectionSelector.H:38
AMREX_FORCE_INLINE IntVect offset(const int face_dir, const int normal)
Definition: ERF_ReadBndryPlanes.cpp:28
Definition: ERF_PlaneAverage.H:14
const amrex::MultiFab & field() const
Definition: ERF_PlaneAverage.H:56
int m_ncell_line
Definition: ERF_PlaneAverage.H:70
AMREX_FORCE_INLINE amrex::Real line_average_interpolated(amrex::Real x, int comp) const
Definition: ERF_PlaneAverage.H:127
int m_precision
Definition: ERF_PlaneAverage.H:72
int ncomp() const
Definition: ERF_PlaneAverage.H:39
amrex::Real m_xlo
Definition: ERF_PlaneAverage.H:67
amrex::IntVect m_ixtype
Definition: ERF_PlaneAverage.H:79
amrex::Vector< amrex::Real > m_line_xcentroid
Definition: ERF_PlaneAverage.H:64
const amrex::Vector< amrex::Real > & line_centroids() const
Definition: ERF_PlaneAverage.H:51
AMREX_FORCE_INLINE void compute_averages(const IndexSelector &idxOp, const amrex::MultiFab &mfab)
void set_precision(int p)
Definition: ERF_PlaneAverage.H:32
const int m_level
Definition: ERF_PlaneAverage.H:73
const amrex::MultiFab * m_field
Definition: ERF_PlaneAverage.H:75
int level() const
Definition: ERF_PlaneAverage.H:38
~PlaneAverage()=default
int m_ncell_plane
Definition: ERF_PlaneAverage.H:69
amrex::Vector< amrex::Real > m_line_average
Definition: ERF_PlaneAverage.H:62
amrex::Real xlo() const
Definition: ERF_PlaneAverage.H:35
AMREX_FORCE_INLINE void operator()()
Definition: ERF_PlaneAverage.H:161
int ncell_plane() const
Definition: ERF_PlaneAverage.H:40
const int m_axis
Definition: ERF_PlaneAverage.H:77
int ncell_line() const
Definition: ERF_PlaneAverage.H:41
amrex::Real dx() const
Definition: ERF_PlaneAverage.H:34
amrex::Real m_dx
Definition: ERF_PlaneAverage.H:66
PlaneAverage()=delete
amrex::Geometry m_geom
Definition: ERF_PlaneAverage.H:76
int axis() const
Definition: ERF_PlaneAverage.H:37
int m_ncomp
Definition: ERF_PlaneAverage.H:59
const amrex::Vector< amrex::Real > & line_average() const
Definition: ERF_PlaneAverage.H:43
amrex::IntVect m_ghost_to_inc
Definition: ERF_PlaneAverage.H:78
@ ng
Definition: ERF_Morrison.H:48