ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
PlaneAverage Class Reference

#include <ERF_PlaneAverage.H>

Collaboration diagram for PlaneAverage:

Public Member Functions

AMREX_FORCE_INLINE PlaneAverage (const amrex::MultiFab *field_in, amrex::Geometry geom_in, int axis_in, amrex::IntVect n_ghost_to_inc=amrex::IntVect{0})
 
 PlaneAverage ()=delete
 
 ~PlaneAverage ()=default
 
AMREX_FORCE_INLINE void operator() ()
 
AMREX_FORCE_INLINE amrex::Real line_average_interpolated (amrex::Real x, int comp) const
 
void set_precision (int p)
 
amrex::Real dx () const
 
amrex::Real xlo () const
 
int axis () const
 
int level () const
 
int ncomp () const
 
int ncell_plane () const
 
int ncell_line () const
 
const amrex::Vector< amrex::Real > & line_average () const
 
AMREX_FORCE_INLINE void line_average (int comp, amrex::Gpu::HostVector< amrex::Real > &l_vec)
 
const amrex::Vector< amrex::Real > & line_centroids () const
 
const amrex::MultiFab & field () const
 
template<typename IndexSelector >
AMREX_FORCE_INLINE void compute_averages (const IndexSelector &idxOp, const amrex::MultiFab &mfab)
 
template<typename IndexSelector >
void compute_averages (const IndexSelector &idxOp, const amrex::MultiFab &mfab)
 

Protected Attributes

int m_ncomp
 
amrex::Vector< amrex::Realm_line_average
 
amrex::Vector< amrex::Realm_line_xcentroid
 
amrex::Real m_dx
 
amrex::Real m_xlo
 
int m_ncell_plane
 
int m_ncell_line
 
int m_precision = 4
 
const int m_level = 0
 
const amrex::MultiFab * m_field
 
amrex::Geometry m_geom
 
const int m_axis
 
amrex::IntVect m_ghost_to_inc = amrex::IntVect{0}
 
amrex::IntVect m_ixtype
 

Detailed Description

Basic averaging and interpolation operations

Constructor & Destructor Documentation

◆ PlaneAverage() [1/2]

PlaneAverage::PlaneAverage ( const amrex::MultiFab *  field_in,
amrex::Geometry  geom_in,
int  axis_in,
amrex::IntVect  n_ghost_to_inc = amrex::IntVect{0} 
)
explicit
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 
109  // Domain guarantees enough space and we usually average in z-dir (spans domain)
111 
112  // First estimate is with domain (updated with MFiter later)
113  m_ncell_plane = 1;
114  auto period = m_geom.periodicity();
115  for (int i = 0; i < AMREX_SPACEDIM; ++i) {
116  int p_fac = (!period.isPeriodic(i)) ? 1 : 0;
117  if (i != m_axis) m_ncell_plane *= (dom_hi[i] - dom_lo[i] + 1 + p_fac*m_ixtype[i]);
118  }
119 
120  m_line_average.resize(static_cast<size_t>(m_ncell_line) * m_ncomp, 0.0);
122 
123  for (int i = 0; i < m_ncell_line; ++i) {
124  m_line_xcentroid[i] = m_xlo + (i + 0.5) * m_dx;
125  }
126 }
AMREX_ALWAYS_ASSERT(bx.length()[2]==khi+1)
const auto & dom_hi
Definition: ERF_SetupVertDiff.H:2
const auto & dom_lo
Definition: ERF_SetupVertDiff.H:1
int m_ncell_line
Definition: ERF_PlaneAverage.H:70
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::MultiFab * m_field
Definition: ERF_PlaneAverage.H:75
int m_ncell_plane
Definition: ERF_PlaneAverage.H:69
amrex::Vector< amrex::Real > m_line_average
Definition: ERF_PlaneAverage.H:62
const int m_axis
Definition: ERF_PlaneAverage.H:77
amrex::Real m_dx
Definition: ERF_PlaneAverage.H:66
amrex::Geometry m_geom
Definition: ERF_PlaneAverage.H:76
int m_ncomp
Definition: ERF_PlaneAverage.H:59
amrex::IntVect m_ghost_to_inc
Definition: ERF_PlaneAverage.H:78
Here is the call graph for this function:

◆ PlaneAverage() [2/2]

PlaneAverage::PlaneAverage ( )
delete

◆ ~PlaneAverage()

PlaneAverage::~PlaneAverage ( )
default

Member Function Documentation

◆ axis()

int PlaneAverage::axis ( ) const
inline
37 { return m_axis; }

◆ compute_averages() [1/2]

template<typename IndexSelector >
AMREX_FORCE_INLINE void PlaneAverage::compute_averages ( const IndexSelector &  idxOp,
const amrex::MultiFab &  mfab 
)

fill line storage with averages

Referenced by Morrison::Compute_Coefficients(), SAM::Compute_Coefficients(), make_mom_sources(), make_sources(), and operator()().

Here is the caller graph for this function:

◆ compute_averages() [2/2]

template<typename IndexSelector >
void PlaneAverage::compute_averages ( const IndexSelector &  idxOp,
const amrex::MultiFab &  mfab 
)
186 {
187  amrex::Vector<amrex::Real> cnt_h(m_ncell_line,0.);
188  amrex::AsyncArray<amrex::Real> cnt_d(cnt_h.data(), cnt_h.size());
189  amrex::AsyncArray<amrex::Real> lavg(m_line_average.data(), m_line_average.size());
190  amrex::Real* cnt_avg = cnt_d.data();
191  amrex::Real* line_avg = lavg.data();
192  const int ncomp = m_ncomp;
193 
194  amrex::Box domain = amrex::convert(m_geom.Domain(),m_ixtype);
195 
196  amrex::IntVect ng = amrex::IntVect(0);
198  ng[m_axis] = offset;
199 
200  std::unique_ptr<amrex::iMultiFab> mask = OwnerMask(*m_field, m_geom.periodicity(), ng);
201 
202 #ifdef _OPENMP
203 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
204 #endif
205  for (amrex::MFIter mfi(mfab, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) {
206  amrex::Box tbx = mfi.tilebox();
207  if (tbx.smallEnd(m_axis) == domain.smallEnd(m_axis)) tbx.growLo(m_axis,offset);
208  if (tbx.bigEnd (m_axis) == domain.bigEnd (m_axis)) tbx.growHi(m_axis,offset);
209  amrex::Box pbx = PerpendicularBox<IndexSelector>(tbx, amrex::IntVect(0));
210 
211  const amrex::Array4<const amrex::Real>& fab_arr = mfab.const_array(mfi);
212  const amrex::Array4<const int >& mask_arr = mask->const_array(mfi);
213 
214  amrex::ParallelFor(amrex::Gpu::KernelInfo().setReduction(true), pbx, [=]
215  AMREX_GPU_DEVICE( int p_i, int p_j, int p_k,
216  amrex::Gpu::Handler const& handler) noexcept
217  {
218  // Loop over the direction perpendicular to the plane.
219  // This reduces the atomic pressure on the destination arrays.
220  amrex::Box lbx = ParallelBox<IndexSelector>(tbx, amrex::IntVect{p_i, p_j, p_k});
221 
222  for (int k = lbx.smallEnd(2); k <= lbx.bigEnd(2); ++k) {
223  for (int j = lbx.smallEnd(1); j <= lbx.bigEnd(1); ++j) {
224  for (int i = lbx.smallEnd(0); i <= lbx.bigEnd(0); ++i) {
225  int ind = idxOp.getIndx(i, j, k) + offset;
226  // NOTE: This factor is to avoid an if statement that will break
227  // the devicereducesum since all threads won't participate.
228  // This more performant than Gpu::Atomic::Add.
229  amrex::Real fac = (mask_arr(i,j,k)) ? 1.0 : 0.0;
230  amrex::Gpu::deviceReduceSum(&cnt_avg[ind], fac, handler);
231  for (int n = 0; n < ncomp; ++n) {
232  amrex::Gpu::deviceReduceSum(&line_avg[ncomp * ind + n],
233  fab_arr(i, j, k, n) * fac, handler);
234  }
235  }
236  }
237  }
238  });
239  }
240 
241  cnt_d.copyToHost(cnt_h.data(), cnt_h.size());
242  lavg.copyToHost(m_line_average.data(), m_line_average.size());
243  amrex::ParallelDescriptor::ReduceRealSum(cnt_h.data(), cnt_h.size());
244  amrex::ParallelDescriptor::ReduceRealSum(m_line_average.data(), m_line_average.size());
245  for (int ind(0); ind<m_ncell_line; ++ind) {
246  for (int n(0); n<ncomp; ++n) {
247  m_line_average[ncomp*ind + n] /= cnt_h[ind];
248  }
249  }
250  m_ncell_plane = static_cast<int>(cnt_h[0]);
251 }
ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { const auto prob_lo=geomdata.ProbLo();const auto dx=geomdata.CellSize();const Real x=(prob_lo[0]+(i+0.5) *dx[0])/mf_m(i, j, 0);const Real z=z_cc(i, j, k);Real L=std::sqrt(std::pow((x - x_c)/x_r, 2)+std::pow((z - z_c)/z_r, 2));if(L<=1.0) { Real dT=T_pert *(std::cos(PI *L)+1.0)/2.0;Real Tbar_hse=p_hse(i, j, k)/(R_d *r_hse(i, j, k));Real theta_perturbed=(Tbar_hse+dT) *std::pow(p_0/p_hse(i, j, k), rdOcp);Real theta_0=(Tbar_hse) *std::pow(p_0/p_hse(i, j, k), rdOcp);if(const_rho) { state_pert(i, j, k, RhoTheta_comp)=r_hse(i, j, k) *(theta_perturbed - theta_0);} else { state_pert(i, j, k, Rho_comp)=getRhoThetagivenP(p_hse(i, j, k))/theta_perturbed - r_hse(i, j, k);} } })
AMREX_FORCE_INLINE IntVect offset(const int face_dir, const int normal)
Definition: ERF_ReadBndryPlanes.cpp:28
amrex::Real Real
Definition: ERF_ShocInterface.H:19
int ncomp() const
Definition: ERF_PlaneAverage.H:39
@ ng
Definition: ERF_Morrison.H:48
Here is the call graph for this function:

◆ dx()

amrex::Real PlaneAverage::dx ( ) const
inline
34 { return m_dx; }

◆ field()

const amrex::MultiFab& PlaneAverage::field ( ) const
inline
56 { return *m_field; }

Referenced by Morrison::Compute_Coefficients(), SAM::Compute_Coefficients(), make_mom_sources(), and make_sources().

Here is the caller graph for this function:

◆ level()

int PlaneAverage::level ( ) const
inline
38 { return m_level; }
const int m_level
Definition: ERF_PlaneAverage.H:73

◆ line_average() [1/2]

const amrex::Vector<amrex::Real>& PlaneAverage::line_average ( ) const
inline
44  {
45  return m_line_average;
46  }

Referenced by Morrison::Compute_Coefficients(), SAM::Compute_Coefficients(), make_mom_sources(), and make_sources().

Here is the caller graph for this function:

◆ line_average() [2/2]

void PlaneAverage::line_average ( int  comp,
amrex::Gpu::HostVector< amrex::Real > &  l_vec 
)
155 {
156  AMREX_ALWAYS_ASSERT(comp >= 0 && comp < m_ncomp);
157 
158  for (int i = 0; i < m_ncell_line; i++) {
159  l_vec[i] = m_line_average[m_ncomp * i + comp];
160  }
161 }
Here is the call graph for this function:

◆ line_average_interpolated()

amrex::Real PlaneAverage::line_average_interpolated ( amrex::Real  x,
int  comp 
) const

evaluate line average at specific location for any average component

130 {
131  AMREX_ALWAYS_ASSERT(comp >= 0 && comp < m_ncomp);
132 
133  amrex::Real c = 0.0;
134  int ind = 0;
135 
136  if (x > m_xlo + 0.5 * m_dx) {
137  ind = static_cast<int>(floor((x - m_xlo) / m_dx - 0.5));
138  const amrex::Real x1 = m_xlo + (ind + 0.5) * m_dx;
139  c = (x - x1) / m_dx;
140  }
141 
142  if (ind + 1 >= m_ncell_line) {
143  ind = m_ncell_line - 2;
144  c = 1.0;
145  }
146 
147  AMREX_ALWAYS_ASSERT(ind >= 0 && ind + 1 < m_ncell_line);
148 
149  return m_line_average[m_ncomp * ind + comp] * (1.0 - c) +
150  m_line_average[m_ncomp * (ind + 1) + comp] * c;
151 }
Here is the call graph for this function:

◆ line_centroids()

const amrex::Vector<amrex::Real>& PlaneAverage::line_centroids ( ) const
inline
52  {
53  return m_line_xcentroid;
54  }

◆ ncell_line()

int PlaneAverage::ncell_line ( ) const
inline
41 { return m_ncell_line; }

Referenced by Morrison::Compute_Coefficients(), SAM::Compute_Coefficients(), make_mom_sources(), and make_sources().

Here is the caller graph for this function:

◆ ncell_plane()

int PlaneAverage::ncell_plane ( ) const
inline
40 { return m_ncell_plane; }

◆ ncomp()

int PlaneAverage::ncomp ( ) const
inline
39 { return m_ncomp; }

Referenced by compute_averages().

Here is the caller graph for this function:

◆ operator()()

void PlaneAverage::operator() ( )
165 {
166  std::fill(m_line_average.begin(), m_line_average.end(), 0.0);
167  switch (m_axis) {
168  case 0:
170  break;
171  case 1:
173  break;
174  case 2:
176  break;
177  default:
178  amrex::Abort("axis must be equal to 0, 1, or 2");
179  break;
180  }
181 }
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 void compute_averages(const IndexSelector &idxOp, const amrex::MultiFab &mfab)
Here is the call graph for this function:

◆ set_precision()

void PlaneAverage::set_precision ( int  p)
inline

change precision of text file output

32 { m_precision = p; }
Real * p
Definition: ERF_InitCustomPert_SquallLine.H:61
int m_precision
Definition: ERF_PlaneAverage.H:72

◆ xlo()

amrex::Real PlaneAverage::xlo ( ) const
inline
35 { return m_xlo; }

Member Data Documentation

◆ m_axis

const int PlaneAverage::m_axis
protected

◆ m_dx

amrex::Real PlaneAverage::m_dx
protected

line storage for centroids of each cell along a line

Referenced by dx(), line_average_interpolated(), and PlaneAverage().

◆ m_field

const amrex::MultiFab* PlaneAverage::m_field
protected

level for plane averaging for now fixed at level=0

Referenced by compute_averages(), field(), operator()(), and PlaneAverage().

◆ m_geom

amrex::Geometry PlaneAverage::m_geom
protected

Referenced by compute_averages(), and PlaneAverage().

◆ m_ghost_to_inc

amrex::IntVect PlaneAverage::m_ghost_to_inc = amrex::IntVect{0}
protected

Referenced by compute_averages(), and PlaneAverage().

◆ m_ixtype

amrex::IntVect PlaneAverage::m_ixtype
protected

Referenced by compute_averages(), and PlaneAverage().

◆ m_level

const int PlaneAverage::m_level = 0
protected

precision for line plot text file

Referenced by level().

◆ m_line_average

amrex::Vector<amrex::Real> PlaneAverage::m_line_average
protected

number of average components line storage for the average velocity and tracer variables

Referenced by compute_averages(), line_average(), line_average_interpolated(), operator()(), and PlaneAverage().

◆ m_line_xcentroid

amrex::Vector<amrex::Real> PlaneAverage::m_line_xcentroid
protected

Referenced by line_centroids(), and PlaneAverage().

◆ m_ncell_line

int PlaneAverage::m_ncell_line
protected

◆ m_ncell_plane

int PlaneAverage::m_ncell_plane
protected

bottom of domain in axis direction

Referenced by compute_averages(), ncell_plane(), and PlaneAverage().

◆ m_ncomp

int PlaneAverage::m_ncomp
protected

◆ m_precision

int PlaneAverage::m_precision = 4
protected

number of cells along line

Referenced by set_precision().

◆ m_xlo

amrex::Real PlaneAverage::m_xlo
protected

mesh spacing in axis direction

Referenced by line_average_interpolated(), PlaneAverage(), and xlo().


The documentation for this class was generated from the following file: