ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
Public Member Functions | Static Public Member Functions | Protected Attributes | List of all members
MOSTAverage Class Reference

#include <ERF_MOSTAverage.H>

Collaboration diagram for MOSTAverage:
Collaboration graph
[legend]

Public Member Functions

 MOSTAverage (amrex::Vector< amrex::Geometry > geom, const bool &has_zphys, std::string a_pp_prefix, const MeshType &m_mesh_type, const TerrainType &m_terrain_type, const amrex::Vector< const eb_ * > &eb_vec={})
 
 ~MOSTAverage ()
 
 MOSTAverage (MOSTAverage &&) noexcept=default
 
MOSTAverage & operator= (MOSTAverage &&other) noexcept=delete
 
 MOSTAverage (const MOSTAverage &other)=delete
 
MOSTAverage & operator= (const MOSTAverage &other)=delete
 
void make_MOSTAverage_at_level (const int &lev, const amrex::Vector< amrex::MultiFab * > &vars_old, std::unique_ptr< amrex::MultiFab > &Theta_prim, std::unique_ptr< amrex::MultiFab > &Qv_prim, std::unique_ptr< amrex::MultiFab > &Qr_prim, std::unique_ptr< amrex::MultiFab > &z_phys_nd)
 
void update_field_ptrs (const int &lev, amrex::Vector< amrex::Vector< amrex::MultiFab >> &vars_old, amrex::Vector< std::unique_ptr< amrex::MultiFab >> &Theta_prim, amrex::Vector< std::unique_ptr< amrex::MultiFab >> &Qv_prim, amrex::Vector< std::unique_ptr< amrex::MultiFab >> &Qr_prim)
 
void set_rotated_fields (const int &lev)
 
void set_plane_normalization (const int &lev)
 
void set_eb_normalization (const int &lev)
 
void set_region_normalization (const int &)
 
void set_k_indices_N (const int &lev)
 
void set_k_indices_T (const int &lev)
 
void set_norm_indices_T (const int &lev)
 
void set_z_positions_T (const int &lev)
 
void set_z_positions_EB (const int &lev)
 
void set_norm_positions_T (const int &lev)
 
void compute_averages (const int &lev)
 
void compute_plane_averages (const int &lev)
 
void compute_region_averages (const int &lev)
 
void compute_eb_averages (const int &lev)
 
void write_k_indices (const int &lev)
 
void write_norm_indices (const int &lev)
 
void write_xz_positions (const int &lev, const int &j)
 
void write_averages (const int &lev)
 
const amrex::MultiFab * get_average (const int &lev, const int &comp) const
 
amrex::MultiFab * get_zref (const int &lev) const
 
const amrex::iMultiFab * get_k_indices (const int &lev) const
 

Static Public Member Functions

AMREX_GPU_HOST_DEVICE static AMREX_INLINE void trilinear_interp_T (const amrex::Real &xp, const amrex::Real &yp, const amrex::Real &zp, amrex::Real *interp_vals, amrex::Array4< amrex::Real const > const &interp_array, amrex::Array4< amrex::Real const > const &z_arr, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &plo, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &dxi, const int interp_comp)
 

Protected Attributes

const amrex::Vector< amrex::Geometry > m_geom
 
amrex::Vector< amrex::Vector< amrex::MultiFab * > > m_fields
 
amrex::Vector< amrex::MultiFab * > m_z_phys_nd
 
std::string m_pp_prefix
 
MeshType m_mesh_type
 
TerrainType m_terrain_type
 
int m_nvar {6}
 
int m_navg {6}
 
int m_maxlev {0}
 
int m_policy {0}
 
bool m_rotate {false}
 
amrex::Vector< std::unique_ptr< amrex::MultiFab > > m_zref
 
amrex::Vector< std::unique_ptr< amrex::MultiFab > > m_x_pos
 
amrex::Vector< std::unique_ptr< amrex::MultiFab > > m_y_pos
 
amrex::Vector< std::unique_ptr< amrex::MultiFab > > m_z_pos
 
amrex::Vector< std::unique_ptr< amrex::iMultiFab > > m_i_indx
 
amrex::Vector< std::unique_ptr< amrex::iMultiFab > > m_j_indx
 
amrex::Vector< std::unique_ptr< amrex::iMultiFab > > m_k_indx
 
amrex::Vector< amrex::Vector< std::unique_ptr< amrex::MultiFab > > > m_averages
 
amrex::Vector< amrex::Vector< std::unique_ptr< amrex::MultiFab > > > m_rot_fields
 
amrex::Vector< amrex::Vector< int > > m_ncell_plane
 
amrex::Vector< amrex::Vector< amrex::Real > > m_plane_average
 
int m_radius {0}
 
int m_ncell_region {1}
 
amrex::Vector< int > m_k_in
 
bool m_interp {false}
 
bool m_norm_vec {false}
 
amrex::Vector< const eb_ * > m_eb_vec
 
amrex::Vector< amrex::Vector< amrex::Real > > m_total_bndry_area
 
bool m_t_avg {false}
 
amrex::Vector< int > m_t_init
 
amrex::Real m_time_window {amrex::Real(1.0e-16)}
 
amrex::Real m_fact_new
 
amrex::Real m_fact_old
 
bool include_subgrid_vel = false
 
amrex::Vector< amrex::Real > m_Vsg
 
const amrex::Real zref_default = amrex::Real(10.0)
 

Constructor & Destructor Documentation

◆ MOSTAverage() [1/3]

MOSTAverage::MOSTAverage ( amrex::Vector< amrex::Geometry >  geom,
const bool &  has_zphys,
std::string  a_pp_prefix,
const MeshType &  m_mesh_type,
const TerrainType &  m_terrain_type,
const amrex::Vector< const eb_ * > &  eb_vec = {} 
)
explicit

◆ ~MOSTAverage()

MOSTAverage::~MOSTAverage ( )
inline
26  {}

◆ MOSTAverage() [2/3]

MOSTAverage::MOSTAverage ( MOSTAverage &&  )
defaultnoexcept

◆ MOSTAverage() [3/3]

MOSTAverage::MOSTAverage ( const MOSTAverage &  other)
delete

Member Function Documentation

◆ compute_averages()

void MOSTAverage::compute_averages ( const int &  lev)

Function to call the type of average computation.

Parameters
[in]levCurrent level
901 {
902  if (m_rotate) set_rotated_fields(lev);
903 
904  switch(m_policy) {
905  case 0: // Standard plane average
906  compute_plane_averages(lev);
907  break;
908  case 1: // Local region/point
909  compute_region_averages(lev);
910  break;
911  case 2: // EB Terrain average
912  compute_eb_averages(lev);
913  break;
914  default:
915  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(false, "Unknown policy for MOSTAverage!");
916  }
917 
918  // We have initialized the averages
919  if (m_t_avg) m_t_init[lev] = 1;
920 }
MOSTAverage::m_t_avg
bool m_t_avg
Definition: ERF_MOSTAverage.H:250
MOSTAverage::compute_plane_averages
void compute_plane_averages(const int &lev)
Definition: ERF_MOSTAverage.cpp:929
MOSTAverage::m_policy
int m_policy
Definition: ERF_MOSTAverage.H:215
MOSTAverage::compute_region_averages
void compute_region_averages(const int &lev)
Definition: ERF_MOSTAverage.cpp:1228
MOSTAverage::m_t_init
amrex::Vector< int > m_t_init
Definition: ERF_MOSTAverage.H:251
MOSTAverage::set_rotated_fields
void set_rotated_fields(const int &lev)
Definition: ERF_MOSTAverage.cpp:298
MOSTAverage::m_rotate
bool m_rotate
Definition: ERF_MOSTAverage.H:216
MOSTAverage::compute_eb_averages
void compute_eb_averages(const int &lev)
Definition: ERF_MOSTAverage.cpp:1605
Here is the call graph for this function:

◆ compute_eb_averages()

void MOSTAverage::compute_eb_averages ( const int &  lev)

Function to compute average over an EB surface.

Parameters
[in]levCurrent level
1606 {
1607  AMREX_ALWAYS_ASSERT(m_eb_vec[lev] != nullptr);
1608 
1609  // Peel back the level
1610  auto& fields = m_fields[lev];
1611  auto& averages = m_averages[lev];
1612  auto& plane_average = m_plane_average[lev];
1613 
1614  // Get EB data
1615  const auto& cc_flags = m_eb_vec[lev]->get_const_factory()->getMultiEBCellFlagFab();
1616  auto cc_afrac = m_eb_vec[lev]->get_const_factory()->getAreaFrac();
1617  const auto& cc_bnorm = m_eb_vec[lev]->get_const_factory()->getBndryNormal();
1618  const auto& u_vfrac = m_eb_vec[lev]->get_u_const_factory()->getVolFrac();
1619  const auto& v_vfrac = m_eb_vec[lev]->get_v_const_factory()->getVolFrac();
1620  const auto& w_vfrac = m_eb_vec[lev]->get_w_const_factory()->getVolFrac();
1621 
1622  // Get geometry for cell sizes
1623  auto const& dx_arr = m_geom[lev].CellSizeArray();
1624  Real dx = dx_arr[0];
1625  Real dy = dx_arr[1];
1626  Real dz = dx_arr[2];
1627 
1628  // Set factors for time averaging
1629  Real d_fact_new, d_fact_old;
1630  if (m_t_avg && m_t_init[lev]) {
1631  d_fact_new = m_fact_new;
1632  d_fact_old = m_fact_old;
1633  } else {
1634  d_fact_new = one;
1635  d_fact_old = zero;
1636  }
1637 
1638  // GPU array to accumulate averages into
1639  Gpu::DeviceVector<Real> pavg(plane_average.size(), zero);
1640  Real* plane_avg = pavg.data();
1641 
1642  // Vectors for normalization and buffer storage
1643  Vector<Real> denom(plane_average.size(),zero);
1644  Vector<Real> val_old(plane_average.size(),zero);
1645 
1646  //
1647  //----------------------------------------------------------
1648  // Averages for U, V, and tangential velocity (in local coordinate)
1649  //----------------------------------------------------------
1650  //
1651  {
1652  denom[0] = one / m_total_bndry_area[lev][0];
1653  val_old[0] = plane_average[0]*d_fact_old;
1654  denom[1] = one / m_total_bndry_area[lev][1];
1655  val_old[1] = plane_average[1]*d_fact_old;
1656 
1657  int iavg = m_navg - 1; // Tangential velocity magnitude
1658  denom[iavg] = one / m_total_bndry_area[lev][iavg];
1659  val_old[iavg] = plane_average[iavg]*d_fact_old;
1660 
1661  const Real Vsg = m_Vsg[lev]; // Subgrid scale velocity
1662 
1663 #ifdef _OPENMP
1664 #pragma omp parallel if (Gpu::notInLaunchRegion())
1665 #endif
1666  for (MFIter mfi(*fields[2], TileNoZ()); mfi.isValid(); ++mfi) {
1667  const auto& flag = cc_flags[mfi];
1668 
1669  // Skip boxes that are not singlevalued (MultiCutFab only has data for singlevalued boxes)
1670  if (flag.getType() != FabType::singlevalued) continue;
1671 
1672  Box bx = mfi.tilebox(); // Full 3D box
1673 
1674  // Get EB arrays
1675  auto const flag_arr = flag.const_array();
1676  auto const afrac_x = cc_afrac[0]->const_array(mfi);
1677  auto const afrac_y = cc_afrac[1]->const_array(mfi);
1678  auto const afrac_z = cc_afrac[2]->const_array(mfi);
1679  auto const bnorm_arr = cc_bnorm.const_array(mfi);
1680  auto const u_vf_arr = u_vfrac.const_array(mfi);
1681  auto const v_vf_arr = v_vfrac.const_array(mfi);
1682  auto const w_vf_arr = w_vfrac.const_array(mfi);
1683 
1684  // Get velocity arrays
1685  auto const u_arr = fields[0]->const_array(mfi);
1686  auto const v_arr = fields[1]->const_array(mfi);
1687  auto const w_arr = fields[5]->const_array(mfi);
1688 
1689  ParallelFor(Gpu::KernelInfo().setReduction(true), bx, [=]
1690  AMREX_GPU_DEVICE(int i, int j, int k, Gpu::Handler const& handler) noexcept
1691  {
1692  // Area-weighted averaging over cut cells at any k
1693  if (flag_arr(i,j,k).isSingleValued()) {
1694  // Compute area from face-centered area fractions
1695  Real axm = afrac_x(i ,j ,k );
1696  Real axp = afrac_x(i+1,j ,k );
1697  Real aym = afrac_y(i ,j ,k );
1698  Real ayp = afrac_y(i ,j+1,k );
1699  Real azm = afrac_z(i ,j ,k );
1700  Real azp = afrac_z(i ,j ,k+1);
1701 
1702  Real adx = (axm - axp) * dy * dz;
1703  Real ady = (aym - ayp) * dx * dz;
1704  Real adz = (azm - azp) * dx * dy;
1705 
1706  Real area = std::sqrt(adx*adx + ady*ady + adz*adz);
1707 
1708  // Volume-weighted interpolation of velocities to cell center
1709  Real vf_u_lo = u_vf_arr(i,j,k);
1710  Real vf_u_hi = u_vf_arr(i+1,j,k);
1711  Real sum_vf_u = vf_u_lo + vf_u_hi;
1712  Real u_cc = (sum_vf_u > zero) ? (u_arr(i,j,k)*vf_u_lo + u_arr(i+1,j,k)*vf_u_hi) / sum_vf_u : zero;
1713 
1714  Real vf_v_lo = v_vf_arr(i,j,k);
1715  Real vf_v_hi = v_vf_arr(i,j+1,k);
1716  Real sum_vf_v = vf_v_lo + vf_v_hi;
1717  Real v_cc = (sum_vf_v > zero) ? (v_arr(i,j,k)*vf_v_lo + v_arr(i,j+1,k)*vf_v_hi) / sum_vf_v : zero;
1718 
1719  Real vf_w_lo = w_vf_arr(i,j,k);
1720  Real vf_w_hi = w_vf_arr(i,j,k+1);
1721  Real sum_vf_w = vf_w_lo + vf_w_hi;
1722  Real w_cc = (sum_vf_w > zero) ? (w_arr(i,j,k)*vf_w_lo + w_arr(i,j,k+1)*vf_w_hi) / sum_vf_w : zero;
1723 
1724  // Get normal vector components
1725  Real nx = bnorm_arr(i,j,k,0);
1726  Real ny = bnorm_arr(i,j,k,1);
1727  Real nz = bnorm_arr(i,j,k,2);
1728 
1729  // Compute tangential velocity components
1730  Real v_dot_n = u_cc*nx + v_cc*ny + w_cc*nz;
1731  Real u_tangent = u_cc - v_dot_n * nx;
1732  Real v_tangent = v_cc - v_dot_n * ny;
1733  Real mag = std::sqrt(u_tangent*u_tangent + v_tangent*v_tangent + Vsg*Vsg);
1734 
1735  // Area-weighted sum
1736  Real val_u = u_tangent * area;
1737  Real val_v = v_tangent * area;
1738  Real val_mag = mag * area;
1739 
1740  Gpu::deviceReduceSum(&plane_avg[0], val_u, handler);
1741  Gpu::deviceReduceSum(&plane_avg[1], val_v, handler);
1742  Gpu::deviceReduceSum(&plane_avg[iavg], val_mag, handler);
1743  }
1744  });
1745  }
1746  }
1747 
1748  //
1749  //----------------------------------------------------------
1750  // Averages for T,Qv (cell-centered scalars)
1751  //----------------------------------------------------------
1752  //
1753  for (int imf(2); imf < 4; ++imf) {
1754 
1755  // Continue if no valid Qv pointer
1756  if (!fields[imf]) continue;
1757 
1758  denom[imf] = one / m_total_bndry_area[lev][imf];
1759  val_old[imf] = plane_average[imf]*d_fact_old;
1760 
1761 #ifdef _OPENMP
1762 #pragma omp parallel if (Gpu::notInLaunchRegion())
1763 #endif
1764  for (MFIter mfi(*fields[imf], TileNoZ()); mfi.isValid(); ++mfi) {
1765  const auto& flag = cc_flags[mfi];
1766 
1767  // Skip boxes that are not singlevalued (MultiCutFab only has data for singlevalued boxes)
1768  if (flag.getType() != FabType::singlevalued) continue;
1769 
1770  Box bx = mfi.tilebox(); // Full 3D box
1771  auto const flag_arr = flag.const_array();
1772  auto const afrac_x = cc_afrac[0]->const_array(mfi);
1773  auto const afrac_y = cc_afrac[1]->const_array(mfi);
1774  auto const afrac_z = cc_afrac[2]->const_array(mfi);
1775  auto const mf_arr = fields[imf]->const_array(mfi);
1776 
1777  ParallelFor(Gpu::KernelInfo().setReduction(true), bx, [=]
1778  AMREX_GPU_DEVICE(int i, int j, int k, Gpu::Handler const& handler) noexcept
1779  {
1780  // Area-weighted averaging over cut cells at any k
1781  if (flag_arr(i,j,k).isSingleValued()) {
1782  // Compute area from face-centered area fractions
1783  Real axm = afrac_x(i ,j ,k );
1784  Real axp = afrac_x(i+1,j ,k );
1785  Real aym = afrac_y(i ,j ,k );
1786  Real ayp = afrac_y(i ,j+1,k );
1787  Real azm = afrac_z(i ,j ,k );
1788  Real azp = afrac_z(i ,j ,k+1);
1789 
1790  Real adx = (axm - axp) * dy * dz;
1791  Real ady = (aym - ayp) * dx * dz;
1792  Real adz = (azm - azp) * dx * dy;
1793 
1794  Real area = std::sqrt(adx*adx + ady*ady + adz*adz);
1795 
1796  Real val = mf_arr(i,j,k) * area;
1797  Gpu::deviceReduceSum(&plane_avg[imf], val, handler);
1798  }
1799  });
1800  }
1801  }
1802 
1803  //
1804  //------------------------------------------------------------------------
1805  // Averages for virtual potential temperature
1806  //------------------------------------------------------------------------
1807  //
1808  if (fields[3]) // We have water vapor
1809  {
1810  int iavg = 4;
1811  denom[iavg] = one / m_total_bndry_area[lev][iavg];
1812  val_old[iavg] = plane_average[iavg]*d_fact_old;
1813 
1814 #ifdef _OPENMP
1815 #pragma omp parallel if (Gpu::notInLaunchRegion())
1816 #endif
1817  for (MFIter mfi(*fields[3], TileNoZ()); mfi.isValid(); ++mfi)
1818  {
1819  const auto& flag = cc_flags[mfi];
1820 
1821  // Skip boxes that are not singlevalued (MultiCutFab only has data for singlevalued boxes)
1822  if (flag.getType() != FabType::singlevalued) continue;
1823 
1824  Box bx = mfi.tilebox(); // Full 3D box
1825  auto const flag_arr = flag.const_array();
1826  auto const afrac_x = cc_afrac[0]->const_array(mfi);
1827  auto const afrac_y = cc_afrac[1]->const_array(mfi);
1828  auto const afrac_z = cc_afrac[2]->const_array(mfi);
1829 
1830  const Array4<Real const> T_mf_arr = fields[2]->const_array(mfi);
1831  const Array4<Real const> qv_mf_arr = fields[3]->const_array(mfi);
1832  const Array4<Real const> qr_mf_arr = (fields[4]) ? fields[4]->const_array(mfi) :
1833  Array4<const Real> {};
1834 
1835  ParallelFor(Gpu::KernelInfo().setReduction(true), bx, [=]
1836  AMREX_GPU_DEVICE(int i, int j, int k, Gpu::Handler const& handler) noexcept
1837  {
1838  if (flag_arr(i,j,k).isSingleValued()) {
1839  // Compute area from face-centered area fractions
1840  Real axm = afrac_x(i ,j ,k );
1841  Real axp = afrac_x(i+1,j ,k );
1842  Real aym = afrac_y(i ,j ,k );
1843  Real ayp = afrac_y(i ,j+1,k );
1844  Real azm = afrac_z(i ,j ,k );
1845  Real azp = afrac_z(i ,j ,k+1);
1846 
1847  Real adx = (axm - axp) * dy * dz;
1848  Real ady = (aym - ayp) * dx * dz;
1849  Real adz = (azm - azp) * dx * dy;
1850 
1851  Real area = std::sqrt(adx*adx + ady*ady + adz*adz);
1852 
1853  Real vfac;
1854  if (qr_mf_arr) {
1855  // We also have liquid water
1856  vfac = one + Real(0.61)*qv_mf_arr(i,j,k) - qr_mf_arr(i,j,k);
1857  } else {
1858  vfac = one + Real(0.61)*qv_mf_arr(i,j,k);
1859  }
1860  const Real val = T_mf_arr(i,j,k) * vfac * area;
1861  Gpu::deviceReduceSum(&plane_avg[iavg], val, handler);
1862  }
1863  });
1864  }
1865  }
1866  else // copy temperature
1867  {
1868  int iavg = m_navg - 2;
1869  denom[iavg] = one / m_total_bndry_area[lev][iavg];
1870  // plane_avg[iavg] = plane_avg[2]
1871  Gpu::copy(Gpu::deviceToDevice, pavg.begin() + 2, pavg.begin() + 3,
1872  pavg.begin() + iavg);
1873  }
1874 
1875  // Copy to host and sum across procs
1876  Gpu::copy(Gpu::deviceToHost, pavg.begin(), pavg.end(), plane_average.begin());
1877  ParallelDescriptor::ReduceRealSum(plane_average.data(), plane_average.size());
1878 
1879  // Normalize by total area and apply time averaging
1880  for (int iavg(0); iavg < m_navg; ++iavg){
1881  plane_average[iavg] *= denom[iavg]*d_fact_new;
1882  plane_average[iavg] += val_old[iavg];
1883  averages[iavg]->setVal(plane_average[iavg]);
1884  }
1885 }
one
constexpr amrex::Real one
Definition: ERF_Constants.H:7
zero
constexpr amrex::Real zero
Definition: ERF_Constants.H:6
vfac
Real vfac
Definition: ERF_InitCustomPertVels_ABL.H:26
dy
const Real dy
Definition: ERF_InitCustomPert_ABL.H:24
dx
const Real dx
Definition: ERF_InitCustomPert_ABL.H:23
AMREX_ALWAYS_ASSERT
AMREX_ALWAYS_ASSERT(bx.length()[2]==khi+1)
ParallelFor
ParallelFor(grown_box, [=] AMREX_GPU_DEVICE(int i, int j, int k) { qrcuten_arr(i, j, k)=Real(0);qscuten_arr(i, j, k)=Real(0);qicuten_arr(i, j, k)=Real(0);})
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:19
TileNoZ
AMREX_FORCE_INLINE amrex::IntVect TileNoZ()
Definition: ERF_TileNoZ.H:11
MOSTAverage::m_navg
int m_navg
Definition: ERF_MOSTAverage.H:213
MOSTAverage::m_averages
amrex::Vector< amrex::Vector< std::unique_ptr< amrex::MultiFab > > > m_averages
Definition: ERF_MOSTAverage.H:224
MOSTAverage::m_total_bndry_area
amrex::Vector< amrex::Vector< amrex::Real > > m_total_bndry_area
Definition: ERF_MOSTAverage.H:246
MOSTAverage::m_Vsg
amrex::Vector< amrex::Real > m_Vsg
Definition: ERF_MOSTAverage.H:258
MOSTAverage::m_plane_average
amrex::Vector< amrex::Vector< amrex::Real > > m_plane_average
Definition: ERF_MOSTAverage.H:230
MOSTAverage::m_fact_new
amrex::Real m_fact_new
Definition: ERF_MOSTAverage.H:253
MOSTAverage::m_fields
amrex::Vector< amrex::Vector< amrex::MultiFab * > > m_fields
Definition: ERF_MOSTAverage.H:204
MOSTAverage::m_eb_vec
amrex::Vector< const eb_ * > m_eb_vec
Definition: ERF_MOSTAverage.H:245
MOSTAverage::m_fact_old
amrex::Real m_fact_old
Definition: ERF_MOSTAverage.H:253
MOSTAverage::m_geom
const amrex::Vector< amrex::Geometry > m_geom
Definition: ERF_MOSTAverage.H:203
WSM6SedCellScratch::dz
@ dz
Definition: ERF_AdvanceWSM6.cpp:104

Referenced by compute_averages().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ compute_plane_averages()

void MOSTAverage::compute_plane_averages ( const int &  lev)

Function to compute average over a plane.

Parameters
[in]levCurrent level
930 {
931  // Peel back the level
932  auto& fields = m_fields[lev];
933  auto& rot_fields = m_rot_fields[lev];
934  auto& averages = m_averages[lev];
935  const auto & geom = m_geom[lev];
936 
937  auto& z_phys = m_z_phys_nd[lev];
938  auto& x_pos = m_x_pos[lev];
939  auto& y_pos = m_y_pos[lev];
940  auto& z_pos = m_z_pos[lev];
941 
942  auto& i_indx = m_i_indx[lev];
943  auto& j_indx = m_j_indx[lev];
944  auto& k_indx = m_k_indx[lev];
945 
946  auto& ncell_plane = m_ncell_plane[lev];
947  auto& plane_average = m_plane_average[lev];
948 
949  // Set factors for time averaging
950  Real d_fact_new, d_fact_old;
951  if (m_t_avg && m_t_init[lev]) {
952  d_fact_new = m_fact_new;
953  d_fact_old = m_fact_old;
954  } else {
955  d_fact_new = one;
956  d_fact_old = zero;
957  }
958 
959  int klo = m_geom[lev].Domain().smallEnd(2);
960 
961  // GPU array to accumulate averages into
962  Gpu::DeviceVector<Real> pavg(plane_average.size(), zero);
963  Real* plane_avg = pavg.data();
964 
965  // Vectors for normalization and buffer storage
966  Vector<Real> denom(plane_average.size(),zero);
967  Vector<Real> val_old(plane_average.size(),zero);
968 
969  //
970  //----------------------------------------------------------
971  // Averages over all the fields
972  //----------------------------------------------------------
973  //
974  Box domain = geom.Domain();
975 
976  Array<int,AMREX_SPACEDIM> is_per = {0,0,0};
977  for (int idim(0); idim < AMREX_SPACEDIM-1; ++idim) {
978  if (geom.isPeriodic(idim)) is_per[idim] = 1;
979  }
980 
981  // Averages for U,V,T,Qv (not Qc or W)
982  for (int imf(0); imf < 4; ++imf) {
983 
984  // Continue if no valid Qv pointer
985  if (!fields[imf]) continue;
986 
987  denom[imf] = one / (Real)ncell_plane[imf];
988  val_old[imf] = plane_average[imf]*d_fact_old;
989 
990 #ifdef _OPENMP
991 #pragma omp parallel if (Gpu::notInLaunchRegion())
992 #endif
993  for (MFIter mfi(*fields[imf], TileNoZ()); mfi.isValid(); ++mfi) {
994  Box vbx = mfi.validbox(); // This is the grid (not tile)
995  Box pbx = mfi.tilebox(); // This is the tile (not grid)
996 
997  if (pbx.smallEnd(2) != klo) { continue; }
998 
999  // Make planar since mfiter is over fields
1000  pbx.makeSlab(2,klo);
1001 
1002  // Avoid double counting nodal data by changing the high end when we are
1003  // at the high side of the grid (not just of the tile)
1004  IndexType ixt = averages[imf]->boxArray().ixType();
1005  for (int idim(0); idim < AMREX_SPACEDIM-1; ++idim) {
1006  if ( ixt.nodeCentered(idim) && (pbx.bigEnd(idim) == vbx.bigEnd(idim)) ) {
1007  int dom_hi = domain.bigEnd(idim)+1;
1008  if (pbx.bigEnd(idim) < dom_hi || is_per[idim]) {
1009  pbx.growHi(idim,-1);
1010  }
1011  }
1012  }
1013 
1014  auto mf_arr = (m_rotate) ? rot_fields[imf]->const_array(mfi) :
1015  fields[imf]->const_array(mfi);
1016 
1017  if (m_interp) {
1018  const auto plo = geom.ProbLoArray();
1019  const auto dxInv = geom.InvCellSizeArray();
1020  const auto z_phys_arr = z_phys->const_array(mfi);
1021  auto x_pos_arr = x_pos->array(mfi);
1022  auto y_pos_arr = y_pos->array(mfi);
1023  auto z_pos_arr = z_pos->array(mfi);
1024  ParallelFor(Gpu::KernelInfo().setReduction(true), pbx, [=]
1025  AMREX_GPU_DEVICE(int i, int j, int , Gpu::Handler const& handler) noexcept
1026  {
1027  Real interp{0};
1028  trilinear_interp_T(x_pos_arr(i,j,0), y_pos_arr(i,j,0), z_pos_arr(i,j,0),
1029  &interp, mf_arr, z_phys_arr, plo, dxInv, 1);
1030  Real val = interp;
1031  Gpu::deviceReduceSum(&plane_avg[imf], val, handler);
1032  });
1033  } else {
1034  auto k_arr = k_indx->const_array(mfi);
1035  auto j_arr = j_indx ? j_indx->const_array(mfi) : Array4<const int> {};
1036  auto i_arr = i_indx ? i_indx->const_array(mfi) : Array4<const int> {};
1037  ParallelFor(Gpu::KernelInfo().setReduction(true), pbx, [=]
1038  AMREX_GPU_DEVICE(int i, int j, int , Gpu::Handler const& handler) noexcept
1039  {
1040  int mk = k_arr(i,j,0);
1041  int mj = j_arr ? j_arr(i,j,0) : j;
1042  int mi = i_arr ? i_arr(i,j,0) : i;
1043  Real val = mf_arr(mi,mj,mk);
1044  Gpu::deviceReduceSum(&plane_avg[imf], val, handler);
1045  });
1046  }
1047  }
1048  }
1049 
1050  //
1051  //------------------------------------------------------------------------
1052  // Averages for virtual potential temperature
1053  // (This is cell-centered so we don't need to worry about double-counting)
1054  //------------------------------------------------------------------------
1055  //
1056  if (fields[3]) // We have water vapor
1057  {
1058  int iavg = 4;
1059  denom[iavg] = one / (Real)ncell_plane[iavg];
1060  val_old[iavg] = plane_average[iavg]*d_fact_old;
1061 
1062 #ifdef _OPENMP
1063 #pragma omp parallel if (Gpu::notInLaunchRegion())
1064 #endif
1065  for (MFIter mfi(*fields[3], TileNoZ()); mfi.isValid(); ++mfi)
1066  {
1067  Box pbx = mfi.tilebox();
1068 
1069  if (pbx.smallEnd(2) != klo) { continue; }
1070 
1071  pbx.makeSlab(2,klo);
1072 
1073  const Array4<Real const>& T_mf_arr = fields[2]->const_array(mfi);
1074  const Array4<Real const>& qv_mf_arr = fields[3]->const_array(mfi);
1075  const Array4<Real const>& qr_mf_arr = (fields[4]) ? fields[4]->const_array(mfi) :
1076  Array4<const Real> {};
1077 
1078  if (m_interp) {
1079  const auto plo = m_geom[lev].ProbLoArray();
1080  const auto dxInv = m_geom[lev].InvCellSizeArray();
1081  const auto z_phys_arr = z_phys->const_array(mfi);
1082  auto x_pos_arr = x_pos->array(mfi);
1083  auto y_pos_arr = y_pos->array(mfi);
1084  auto z_pos_arr = z_pos->array(mfi);
1085  ParallelFor(Gpu::KernelInfo().setReduction(true), pbx, [=]
1086  AMREX_GPU_DEVICE(int i, int j, int , Gpu::Handler const& handler) noexcept
1087  {
1088  Real T_interp{0};
1089  Real qv_interp{0};
1090  trilinear_interp_T(x_pos_arr(i,j,0), y_pos_arr(i,j,0), z_pos_arr(i,j,0),
1091  &T_interp, T_mf_arr, z_phys_arr, plo, dxInv, 1);
1092  trilinear_interp_T(x_pos_arr(i,j,0), y_pos_arr(i,j,0), z_pos_arr(i,j,0),
1093  &qv_interp, qv_mf_arr, z_phys_arr, plo, dxInv, 1);
1094  Real vfac;
1095  if (qr_mf_arr) {
1096  // We also have liquid water
1097  Real qr_interp{0};
1098  trilinear_interp_T(x_pos_arr(i,j,0), y_pos_arr(i,j,0), z_pos_arr(i,j,0),
1099  &qr_interp, qr_mf_arr, z_phys_arr, plo, dxInv, 1);
1100  vfac = one + Real(0.61)*qv_interp - qr_interp;
1101  } else {
1102  vfac = one + Real(0.61)*qv_interp;
1103  }
1104  const Real val = T_interp * vfac;
1105  Gpu::deviceReduceSum(&plane_avg[iavg], val, handler);
1106  });
1107  } else {
1108  auto k_arr = k_indx->const_array(mfi);
1109  auto j_arr = j_indx ? j_indx->const_array(mfi) : Array4<const int> {};
1110  auto i_arr = i_indx ? i_indx->const_array(mfi) : Array4<const int> {};
1111  ParallelFor(Gpu::KernelInfo().setReduction(true), pbx, [=]
1112  AMREX_GPU_DEVICE(int i, int j, int , Gpu::Handler const& handler) noexcept
1113  {
1114  int mk = k_arr(i,j,0);
1115  int mj = j_arr ? j_arr(i,j,0) : j;
1116  int mi = i_arr ? i_arr(i,j,0) : i;
1117  Real vfac;
1118  if (qr_mf_arr) {
1119  // We also have liquid water
1120  vfac = one + Real(0.61)*qv_mf_arr(mi,mj,mk) - qr_mf_arr(mi,mj,mk);
1121  } else {
1122  vfac = one + Real(0.61)*qv_mf_arr(mi,mj,mk);
1123  }
1124  const Real val = T_mf_arr(mi,mj,mk) * vfac;
1125  Gpu::deviceReduceSum(&plane_avg[iavg], val, handler);
1126  });
1127  }
1128  }
1129  }
1130  else // copy temperature
1131  {
1132  int iavg = m_navg - 2;
1133  denom[iavg] = one / (Real)ncell_plane[iavg];
1134  // plane_avg[iavg] = plane_avg[2]
1135  Gpu::copy(Gpu::deviceToDevice, pavg.begin() + 2, pavg.begin() + 3,
1136  pavg.begin() + iavg);
1137  }
1138 
1139  //
1140  //------------------------------------------------------------------------
1141  // Averages for the tangential velocity magnitude
1142  // (This is cell-centered so we don't need to worry about double-counting)
1143  //------------------------------------------------------------------------
1144  //
1145  {
1146  int imf_cc = 2;
1147  int imf = 0;
1148  int iavg = m_navg - 1;
1149  denom[iavg] = one / (Real)ncell_plane[iavg];
1150  val_old[iavg] = plane_average[iavg]*d_fact_old;
1151 
1152  const Real Vsg = m_Vsg[lev];
1153 
1154 #ifdef _OPENMP
1155 #pragma omp parallel if (Gpu::notInLaunchRegion())
1156 #endif
1157  for (MFIter mfi(*fields[imf_cc], TileNoZ()); mfi.isValid(); ++mfi)
1158  {
1159  Box pbx = mfi.tilebox();
1160 
1161  if (pbx.smallEnd(2) != klo) { continue; }
1162 
1163  pbx.makeSlab(2,klo);
1164 
1165  // Last element is Umag and always cell centered
1166  auto u_mf_arr = (m_rotate) ? rot_fields[imf ]->const_array(mfi) :
1167  fields[imf ]->const_array(mfi);
1168  auto v_mf_arr = (m_rotate) ? rot_fields[imf+1]->const_array(mfi) :
1169  fields[imf+1]->const_array(mfi);
1170 
1171  if (m_interp) {
1172  const auto plo = m_geom[lev].ProbLoArray();
1173  const auto dxInv = m_geom[lev].InvCellSizeArray();
1174  const auto z_phys_arr = z_phys->const_array(mfi);
1175  auto x_pos_arr = x_pos->array(mfi);
1176  auto y_pos_arr = y_pos->array(mfi);
1177  auto z_pos_arr = z_pos->array(mfi);
1178  ParallelFor(Gpu::KernelInfo().setReduction(true), pbx, [=]
1179  AMREX_GPU_DEVICE(int i, int j, int , Gpu::Handler const& handler) noexcept
1180  {
1181  Real u_interp{0};
1182  Real v_interp{0};
1183  trilinear_interp_T(x_pos_arr(i,j,0), y_pos_arr(i,j,0), z_pos_arr(i,j,0),
1184  &u_interp, u_mf_arr, z_phys_arr, plo, dxInv, 1);
1185  trilinear_interp_T(x_pos_arr(i,j,0), y_pos_arr(i,j,0), z_pos_arr(i,j,0),
1186  &v_interp, v_mf_arr, z_phys_arr, plo, dxInv, 1);
1187  const Real val = std::sqrt(u_interp*u_interp + v_interp*v_interp + Vsg*Vsg);
1188  Gpu::deviceReduceSum(&plane_avg[iavg], val, handler);
1189  });
1190  } else {
1191  auto k_arr = k_indx->const_array(mfi);
1192  auto j_arr = j_indx ? j_indx->const_array(mfi) : Array4<const int> {};
1193  auto i_arr = i_indx ? i_indx->const_array(mfi) : Array4<const int> {};
1194  ParallelFor(Gpu::KernelInfo().setReduction(true), pbx, [=]
1195  AMREX_GPU_DEVICE(int i, int j, int , Gpu::Handler const& handler) noexcept
1196  {
1197  int mk = k_arr(i,j,0);
1198  int mj = j_arr ? j_arr(i,j,0) : j;
1199  int mi = i_arr ? i_arr(i,j,0) : i;
1200  const Real u_val = myhalf * (u_mf_arr(mi,mj,mk) + u_mf_arr(mi+1,mj ,mk));
1201  const Real v_val = myhalf * (v_mf_arr(mi,mj,mk) + v_mf_arr(mi ,mj+1,mk));
1202  const Real val = std::sqrt(u_val*u_val + v_val*v_val + Vsg*Vsg);
1203  Gpu::deviceReduceSum(&plane_avg[iavg], val, handler);
1204  });
1205  }
1206  }
1207  }
1208 
1209  // Copy to host and sum across procs
1210  Gpu::copy(Gpu::deviceToHost, pavg.begin(), pavg.end(), plane_average.begin());
1211  ParallelDescriptor::ReduceRealSum(plane_average.data(), static_cast<int>(plane_average.size()));
1212 
1213  // No spatial variation with plane averages
1214  for (int iavg(0); iavg < m_navg; ++iavg){
1215  plane_average[iavg] *= denom[iavg]*d_fact_new;
1216  plane_average[iavg] += val_old[iavg];
1217  averages[iavg]->setVal(plane_average[iavg]);
1218  }
1219 }
myhalf
constexpr amrex::Real myhalf
Definition: ERF_Constants.H:11
dxInv
amrex::GpuArray< Real, AMREX_SPACEDIM > dxInv
Definition: ERF_InitCustomPertVels_ParticleTests.H:17
dom_hi
const auto & dom_hi
Definition: ERF_SetupVertDiff.H:2
MOSTAverage::m_y_pos
amrex::Vector< std::unique_ptr< amrex::MultiFab > > m_y_pos
Definition: ERF_MOSTAverage.H:219
MOSTAverage::m_i_indx
amrex::Vector< std::unique_ptr< amrex::iMultiFab > > m_i_indx
Definition: ERF_MOSTAverage.H:221
MOSTAverage::m_z_phys_nd
amrex::Vector< amrex::MultiFab * > m_z_phys_nd
Definition: ERF_MOSTAverage.H:205
MOSTAverage::m_x_pos
amrex::Vector< std::unique_ptr< amrex::MultiFab > > m_x_pos
Definition: ERF_MOSTAverage.H:218
MOSTAverage::m_rot_fields
amrex::Vector< amrex::Vector< std::unique_ptr< amrex::MultiFab > > > m_rot_fields
Definition: ERF_MOSTAverage.H:225
MOSTAverage::m_z_pos
amrex::Vector< std::unique_ptr< amrex::MultiFab > > m_z_pos
Definition: ERF_MOSTAverage.H:220
MOSTAverage::m_ncell_plane
amrex::Vector< amrex::Vector< int > > m_ncell_plane
Definition: ERF_MOSTAverage.H:229
MOSTAverage::m_j_indx
amrex::Vector< std::unique_ptr< amrex::iMultiFab > > m_j_indx
Definition: ERF_MOSTAverage.H:222
MOSTAverage::m_interp
bool m_interp
Definition: ERF_MOSTAverage.H:240
MOSTAverage::trilinear_interp_T
AMREX_GPU_HOST_DEVICE static AMREX_INLINE void trilinear_interp_T(const amrex::Real &xp, const amrex::Real &yp, const amrex::Real &zp, amrex::Real *interp_vals, amrex::Array4< amrex::Real const > const &interp_array, amrex::Array4< amrex::Real const > const &z_arr, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &plo, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &dxi, const int interp_comp)
Definition: ERF_MOSTAverage.H:134
MOSTAverage::m_k_indx
amrex::Vector< std::unique_ptr< amrex::iMultiFab > > m_k_indx
Definition: ERF_MOSTAverage.H:223

Referenced by compute_averages().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ compute_region_averages()

void MOSTAverage::compute_region_averages ( const int &  lev)

Function to compute average over local region.

Parameters
[in]levCurrent level
1229 {
1230  // Peel back the level
1231  auto& fields = m_fields[lev];
1232  auto& rot_fields = m_rot_fields[lev];
1233  auto& averages = m_averages[lev];
1234  const auto & geom = m_geom[lev];
1235 
1236  auto& z_phys = m_z_phys_nd[lev];
1237  auto& x_pos = m_x_pos[lev];
1238  auto& y_pos = m_y_pos[lev];
1239  auto& z_pos = m_z_pos[lev];
1240 
1241  auto& i_indx = m_i_indx[lev];
1242  auto& j_indx = m_j_indx[lev];
1243  auto& k_indx = m_k_indx[lev];
1244 
1245  int klo = m_geom[lev].Domain().smallEnd(2);
1246 
1247  // Set factors for time averaging
1248  Real d_fact_new, d_fact_old;
1249  if (m_t_avg && m_t_init[lev]) {
1250  d_fact_new = m_fact_new;
1251  d_fact_old = m_fact_old;
1252  } else {
1253  d_fact_new = one;
1254  d_fact_old = zero;
1255  }
1256 
1257  // Number of cells contained in the local average
1258  const Real denom = one / (Real) m_ncell_region;
1259 
1260  // Capture radius for device
1261  int d_radius = m_radius;
1262 
1263  //
1264  //----------------------------------------------------------
1265  // Averages for U,V,T,Qv
1266  //----------------------------------------------------------
1267  //
1268  for (int imf(0); imf < 4; ++imf) {
1269 
1270  // Continue if no valid Qv pointer
1271  if (!fields[imf]) continue;
1272 
1273 #ifdef _OPENMP
1274 #pragma omp parallel if (Gpu::notInLaunchRegion())
1275 #endif
1276  for (MFIter mfi(*fields[imf], TileNoZ()); mfi.isValid(); ++mfi) {
1277  Box pbx = mfi.tilebox();
1278 
1279  if (pbx.smallEnd(2) != klo) { continue; }
1280 
1281  // Make planar since mfiter is over fields
1282  pbx.makeSlab(2,klo);
1283 
1284  auto mf_arr = (m_rotate) ? rot_fields[imf]->const_array(mfi) :
1285  fields[imf]->const_array(mfi);
1286  auto ma_arr = averages[imf]->array(mfi);
1287 
1288  if (m_interp) {
1289  const auto plo = geom.ProbLoArray();
1290  const auto dx = geom.CellSizeArray();
1291  const auto dxInv = geom.InvCellSizeArray();
1292  const auto z_phys_arr = z_phys->const_array(mfi);
1293  auto x_pos_arr = x_pos->array(mfi);
1294  auto y_pos_arr = y_pos->array(mfi);
1295  auto z_pos_arr = z_pos->array(mfi);
1296  ParallelFor(pbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
1297  {
1298  ma_arr(i,j,0) *= d_fact_old;
1299 
1300  Real met_h_zeta = Compute_h_zeta_AtCellCenter(i,j,k,dxInv,z_phys_arr);
1301  for (int lk(-d_radius); lk <= (d_radius); ++lk) {
1302  for (int lj(-d_radius); lj <= (d_radius); ++lj) {
1303  for (int li(-d_radius); li <= (d_radius); ++li) {
1304  Real interp{0};
1305  Real xp = x_pos_arr(i+li,j+lj,0);
1306  Real yp = y_pos_arr(i+li,j+lj,0);
1307  Real zp = z_pos_arr(i+li,j+lj,0) + met_h_zeta*lk*dx[2];
1308  trilinear_interp_T(xp, yp, zp, &interp, mf_arr, z_phys_arr, plo, dxInv, 1);
1309  Real val = denom * interp * d_fact_new;
1310  ma_arr(i,j,0) += val;
1311  }
1312  }
1313  }
1314  });
1315  } else {
1316  auto k_arr = k_indx->const_array(mfi);
1317  auto j_arr = j_indx ? j_indx->const_array(mfi) : Array4<const int> {};
1318  auto i_arr = i_indx ? i_indx->const_array(mfi) : Array4<const int> {};
1319  ParallelFor(pbx, [=] AMREX_GPU_DEVICE(int i, int j, int ) noexcept
1320  {
1321  ma_arr(i,j,0) *= d_fact_old;
1322 
1323  int mk = k_arr(i,j,0);
1324  int mj = j_arr ? j_arr(i,j,0) : j;
1325  int mi = i_arr ? i_arr(i,j,0) : i;
1326  for (int lk(mk-d_radius); lk <= (mk+d_radius); ++lk) {
1327  for (int lj(mj-d_radius); lj <= (mj+d_radius); ++lj) {
1328  for (int li(mi-d_radius); li <= (mi+d_radius); ++li) {
1329  Real val = denom * mf_arr(li, lj, lk) * d_fact_new;
1330  ma_arr(i,j,0) += val;
1331  }
1332  }
1333  }
1334  });
1335  }
1336  } // MFiter
1337 
1338  // Fill interior ghost cells and any ghost cells outside a periodic domain
1339  //***********************************************************************************
1340  averages[imf]->FillBoundary(geom.periodicity());
1341 
1342  } // imf
1343 
1344  //
1345  //----------------------------------------------------------
1346  // Averages for virtual potential temperature
1347  //----------------------------------------------------------
1348  //
1349  if (fields[3]) // We have water vapor
1350  {
1351  int iavg = 4;
1352 
1353 #ifdef _OPENMP
1354 #pragma omp parallel if (Gpu::notInLaunchRegion())
1355 #endif
1356  for (MFIter mfi(*fields[3], TileNoZ()); mfi.isValid(); ++mfi) {
1357  Box pbx = mfi.tilebox();
1358 
1359  if (pbx.smallEnd(2) != klo) { continue; }
1360 
1361  pbx.makeSlab(2,klo);
1362 
1363  const Array4<Real const>& T_mf_arr = fields[2]->const_array(mfi);
1364  const Array4<Real const>& qv_mf_arr = fields[3]->const_array(mfi);
1365  const Array4<Real const>& qr_mf_arr = (fields[4]) ? fields[4]->const_array(mfi) :
1366  Array4<const Real> {};
1367  auto ma_arr = averages[iavg]->array(mfi);
1368 
1369  if (m_interp) {
1370  const auto plo = geom.ProbLoArray();
1371  const auto dx = geom.CellSizeArray();
1372  const auto dxInv = geom.InvCellSizeArray();
1373  const auto z_phys_arr = z_phys->const_array(mfi);
1374  auto x_pos_arr = x_pos->array(mfi);
1375  auto y_pos_arr = y_pos->array(mfi);
1376  auto z_pos_arr = z_pos->array(mfi);
1377  ParallelFor(pbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
1378  {
1379  ma_arr(i,j,0) *= d_fact_old;
1380 
1381  Real met_h_zeta = Compute_h_zeta_AtCellCenter(i,j,k,dxInv,z_phys_arr);
1382  for (int lk(-d_radius); lk <= (d_radius); ++lk) {
1383  for (int lj(-d_radius); lj <= (d_radius); ++lj) {
1384  for (int li(-d_radius); li <= (d_radius); ++li) {
1385  Real T_interp{0};
1386  Real qv_interp{0};
1387  Real xp = x_pos_arr(i+li,j+lj,0);
1388  Real yp = y_pos_arr(i+li,j+lj,0);
1389  Real zp = z_pos_arr(i+li,j+lj,0) + met_h_zeta*lk*dx[2];
1390  trilinear_interp_T(xp, yp, zp, &T_interp, T_mf_arr, z_phys_arr, plo, dxInv, 1);
1391  trilinear_interp_T(xp, yp, zp, &qv_interp, qv_mf_arr, z_phys_arr, plo, dxInv, 1);
1392  Real vfac;
1393  if (qr_mf_arr) {
1394  // We also have liquid water
1395  Real qr_interp{0};
1396  trilinear_interp_T(x_pos_arr(i,j,0), y_pos_arr(i,j,0), z_pos_arr(i,j,0),
1397  &qr_interp, qr_mf_arr, z_phys_arr, plo, dxInv, 1);
1398  vfac = one + Real(0.61)*qv_interp - qr_interp;
1399  } else {
1400  vfac = one + Real(0.61)*qv_interp;
1401  }
1402  const Real mag = T_interp * vfac;
1403  const Real val = denom * mag * d_fact_new;
1404  ma_arr(i,j,0) += val;
1405  }
1406  }
1407  }
1408  });
1409  } else {
1410  auto k_arr = k_indx->const_array(mfi);
1411  auto j_arr = j_indx ? j_indx->const_array(mfi) : Array4<const int> {};
1412  auto i_arr = i_indx ? i_indx->const_array(mfi) : Array4<const int> {};
1413  ParallelFor(pbx, [=] AMREX_GPU_DEVICE(int i, int j, int ) noexcept
1414  {
1415  ma_arr(i,j,0) *= d_fact_old;
1416 
1417  int mk = k_arr(i,j,0);
1418  int mj = j_arr ? j_arr(i,j,0) : j;
1419  int mi = i_arr ? i_arr(i,j,0) : i;
1420  for (int lk(mk-d_radius); lk <= (mk+d_radius); ++lk) {
1421  for (int lj(mj-d_radius); lj <= (mj+d_radius); ++lj) {
1422  for (int li(mi-d_radius); li <= (mi+d_radius); ++li) {
1423  Real vfac;
1424  if (qr_mf_arr) {
1425  // We also have liquid water
1426  vfac = one + Real(0.61)*qv_mf_arr(li,lj,lk) - qr_mf_arr(li,lj,lk);
1427  } else {
1428  vfac = one + Real(0.61)*qv_mf_arr(li,lj,lk);
1429  }
1430  const Real mag = T_mf_arr(li,lj,lk) * vfac;
1431  const Real val = denom * mag * d_fact_new;
1432  ma_arr(i,j,0) += val;
1433  }
1434  }
1435  }
1436  });
1437  }
1438  } // MFiter
1439 
1440  // Fill interior ghost cells and any ghost cells outside a periodic domain
1441  //***********************************************************************************
1442  averages[iavg]->FillBoundary(geom.periodicity());
1443 
1444  }
1445  else // copy temperature
1446  {
1447  int iavg = m_navg - 2;
1448  IntVect ng = averages[iavg]->nGrowVect();
1449  MultiFab::Copy(*(averages[iavg]),*(averages[2]),0,0,1,ng);
1450  }
1451 
1452  //
1453  //----------------------------------------------------------
1454  // Averages for the tangential velocity magnitude
1455  //----------------------------------------------------------
1456  //
1457  {
1458  int imf_cc = 2;
1459  int imf = 0;
1460  int iavg = m_navg - 1;
1461 
1462  const Real Vsg = m_Vsg[lev];
1463 
1464 #ifdef _OPENMP
1465 #pragma omp parallel if (Gpu::notInLaunchRegion())
1466 #endif
1467  for (MFIter mfi(*fields[imf_cc], TileNoZ()); mfi.isValid(); ++mfi) {
1468  Box pbx = mfi.tilebox();
1469 
1470  if (pbx.smallEnd(2) != klo) { continue; }
1471 
1472  pbx.makeSlab(2,klo);
1473 
1474  auto u_mf_arr = (m_rotate) ? rot_fields[imf ]->const_array(mfi) :
1475  fields[imf ]->const_array(mfi);
1476  auto v_mf_arr = (m_rotate) ? rot_fields[imf+1]->const_array(mfi) :
1477  fields[imf+1]->const_array(mfi);
1478  auto ma_arr = averages[iavg]->array(mfi);
1479 
1480  if (m_interp) {
1481  const auto plo = geom.ProbLoArray();
1482  const auto dx = geom.CellSizeArray();
1483  const auto dxInv = geom.InvCellSizeArray();
1484  const auto z_phys_arr = z_phys->const_array(mfi);
1485  auto x_pos_arr = x_pos->array(mfi);
1486  auto y_pos_arr = y_pos->array(mfi);
1487  auto z_pos_arr = z_pos->array(mfi);
1488  ParallelFor(pbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
1489  {
1490  ma_arr(i,j,0) *= d_fact_old;
1491 
1492  Real met_h_zeta = Compute_h_zeta_AtCellCenter(i,j,k,dxInv,z_phys_arr);
1493  for (int lk(-d_radius); lk <= (d_radius); ++lk) {
1494  for (int lj(-d_radius); lj <= (d_radius); ++lj) {
1495  for (int li(-d_radius); li <= (d_radius); ++li) {
1496  Real u_interp{0};
1497  Real v_interp{0};
1498  Real xp = x_pos_arr(i+li,j+lj,0);
1499  Real yp = y_pos_arr(i+li,j+lj,0);
1500  Real zp = z_pos_arr(i+li,j+lj,0) + met_h_zeta*lk*dx[2];
1501  trilinear_interp_T(xp, yp, zp, &u_interp, u_mf_arr, z_phys_arr, plo, dxInv, 1);
1502  trilinear_interp_T(xp, yp, zp, &v_interp, v_mf_arr, z_phys_arr, plo, dxInv, 1);
1503  const Real mag = std::sqrt(u_interp*u_interp + v_interp*v_interp + Vsg*Vsg);
1504  Real val = denom * mag * d_fact_new;
1505  ma_arr(i,j,0) += val;
1506  }
1507  }
1508  }
1509  });
1510  } else {
1511  auto k_arr = k_indx->const_array(mfi);
1512  auto j_arr = j_indx ? j_indx->const_array(mfi) : Array4<const int> {};
1513  auto i_arr = i_indx ? i_indx->const_array(mfi) : Array4<const int> {};
1514  ParallelFor(pbx, [=] AMREX_GPU_DEVICE(int i, int j, int ) noexcept
1515  {
1516  ma_arr(i,j,0) *= d_fact_old;
1517 
1518  int mk = k_arr(i,j,0);
1519  int mj = j_arr ? j_arr(i,j,0) : j;
1520  int mi = i_arr ? i_arr(i,j,0) : i;
1521  for (int lk(mk-d_radius); lk <= (mk+d_radius); ++lk) {
1522  for (int lj(mj-d_radius); lj <= (mj+d_radius); ++lj) {
1523  for (int li(mi-d_radius); li <= (mi+d_radius); ++li) {
1524  const Real u_val = myhalf * (u_mf_arr(li,lj,lk) + u_mf_arr(li+1,lj ,lk));
1525  const Real v_val = myhalf * (v_mf_arr(li,lj,lk) + v_mf_arr(li ,lj+1,lk));
1526  const Real mag = std::sqrt(u_val*u_val + v_val*v_val + Vsg*Vsg);
1527  Real val = denom * mag * d_fact_new;
1528  ma_arr(i,j,0) += val;
1529  }
1530  }
1531  }
1532  });
1533  }
1534  } // MFiter
1535 
1536  // Fill interior ghost cells and any ghost cells outside a periodic domain
1537  //***********************************************************************************
1538  averages[iavg]->FillBoundary(geom.periodicity());
1539 
1540  }
1541 
1542  // NOTE: Checking periodicity with the geom structure is not
1543  // sufficient at higher levels. The BA may be contained
1544  // within the domain and it's exterior ghost cells filled
1545  // from interpolation; yet the domain BCs are periodic.
1546 
1547  // Need to fill ghost cells outside the domain if not periodic
1548  bool not_per_x = !(geom.periodicity().isPeriodic(0));
1549  bool not_per_y = !(geom.periodicity().isPeriodic(1));
1550  Box cc_bnd_bx = (m_fields[lev][2]->boxArray()).minimalBox();
1551  Box domain = geom.Domain();
1552  if (domain.contains(cc_bnd_bx) || (not_per_x || not_per_y)) {
1553  for (int iavg(0); iavg < m_navg; ++iavg) {
1554  IntVect ng = averages[iavg]->nGrowVect(); ng[2]=0;
1555 
1556  // NOTE: Level 0 spans the whole domain, but finer
1557  // levels do not have such a restriction.
1558  // For now, use the bounding box of the boxArray.
1559 
1560  // NOTE2: The fields and averages have different indexing.
1561  // The averages are: U/V/T/Qv/Tv/Umag
1562  // The fields are: U/V/T/Qv/Qr/W
1563  // We clip iavg at 2 since all the remaining data is CC
1564 
1565  // Bounded box of CC data used for normalization
1566  int imf = min(iavg,2);
1567  Box bnd_bx = (fields[imf]->boxArray()).minimalBox();
1568 #ifdef _OPENMP
1569 #pragma omp parallel if (Gpu::notInLaunchRegion())
1570 #endif
1571  for (MFIter mfi(*fields[imf], TileNoZ()); mfi.isValid(); ++mfi) {
1572  Box gpbx = mfi.growntilebox(ng);
1573 
1574  if (gpbx.smallEnd(2) != klo) { continue; }
1575 
1576  gpbx.makeSlab(2,klo);
1577 
1578  if (bnd_bx.contains(gpbx)) continue;
1579 
1580  auto ma_arr = averages[iavg]->array(mfi);
1581 
1582  int i_lo = bnd_bx.smallEnd(0); int i_hi = bnd_bx.bigEnd(0);
1583  int j_lo = bnd_bx.smallEnd(1); int j_hi = bnd_bx.bigEnd(1);
1584  ParallelFor(gpbx, [=] AMREX_GPU_DEVICE(int i, int j, int ) noexcept
1585  {
1586  int li, lj;
1587  li = i < i_lo ? i_lo : i;
1588  li = li > i_hi ? i_hi : li;
1589  lj = j < j_lo ? j_lo : j;
1590  lj = lj > j_hi ? j_hi : lj;
1591 
1592  ma_arr(i,j,0) = ma_arr(li,lj,0);
1593  });
1594  } // MFiter
1595  } // iavg
1596  } // Not periodic
1597 }
Compute_h_zeta_AtCellCenter
AMREX_FORCE_INLINE AMREX_GPU_DEVICE amrex::Real Compute_h_zeta_AtCellCenter(const int &i, const int &j, const int &k, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &z_nd)
Definition: ERF_TerrainMetrics.H:55
MOSTAverage::m_radius
int m_radius
Definition: ERF_MOSTAverage.H:234
MOSTAverage::m_ncell_region
int m_ncell_region
Definition: ERF_MOSTAverage.H:235
MicVar_Morr::ng
@ ng
Definition: ERF_Morrison.H:48

Referenced by compute_averages().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ get_average()

const amrex::MultiFab* MOSTAverage::get_average ( const int &  lev,
const int &  comp 
) const
inline
113 { return m_averages[lev][comp].get(); }

Referenced by SurfaceLayer::get_mac_avg().

Here is the caller graph for this function:

◆ get_k_indices()

const amrex::iMultiFab* MOSTAverage::get_k_indices ( const int &  lev) const
inline
119 { return m_k_indx[lev].get(); }

◆ get_zref()

amrex::MultiFab* MOSTAverage::get_zref ( const int &  lev) const
inline
116 { return m_zref[lev].get(); }
MOSTAverage::m_zref
amrex::Vector< std::unique_ptr< amrex::MultiFab > > m_zref
Definition: ERF_MOSTAverage.H:217

Referenced by SurfaceLayer::get_zref().

Here is the caller graph for this function:

◆ make_MOSTAverage_at_level()

void MOSTAverage::make_MOSTAverage_at_level ( const int &  lev,
const amrex::Vector< amrex::MultiFab * > &  vars_old,
std::unique_ptr< amrex::MultiFab > &  Theta_prim,
std::unique_ptr< amrex::MultiFab > &  Qv_prim,
std::unique_ptr< amrex::MultiFab > &  Qr_prim,
std::unique_ptr< amrex::MultiFab > &  z_phys_nd 
)
96 {
97  m_fields[lev].resize(m_nvar);
98  m_rot_fields[lev].resize(m_nvar-1);
99  m_averages[lev].resize(m_navg);
100  m_z_phys_nd[lev] = z_phys_nd.get();
101 
102  bool use_terrain_fitted_coords = ( (m_terrain_type == TerrainType::StaticFittedMesh) ||
103  (m_terrain_type == TerrainType::MovingFittedMesh) );
104 
105  bool use_eb = (m_terrain_type == TerrainType::EB);
106 
107  { // Nodal in x
108  auto& mf = *vars_old[Vars::xvel];
109  // Create a 2D ba, dm, & ghost cells
110  const BoxArray& ba = mf.boxArray();
111  BoxList bl2d = ba.boxList();
112  for (auto& b : bl2d) { b.setRange(2,0); }
113  BoxArray ba2d(std::move(bl2d));
114  const DistributionMapping& dm = mf.DistributionMap();
115  const int ncomp = 1;
116  IntVect ng = mf.nGrowVect(); ng[2]=0;
117 
118  m_fields[lev][0] = vars_old[Vars::xvel];
119  m_averages[lev][0] = std::make_unique<MultiFab>(ba2d,dm,ncomp,ng);
120  m_averages[lev][0]->setVal(Real(1e34));
121  if (m_rotate) {
122  m_rot_fields[lev][0] = std::make_unique<MultiFab>(ba,dm,ncomp,ng);
123  MultiFab::Copy(*m_rot_fields[lev][0],mf,0,0,1,ng);
124  } else {
125  m_rot_fields[lev][0] = nullptr;
126  }
127  }
128  { // Nodal in y
129  auto& mf = *vars_old[Vars::yvel];
130  // Create a 2D ba, dm, & ghost cells
131  const BoxArray& ba = mf.boxArray();
132  BoxList bl2d = ba.boxList();
133  for (auto& b : bl2d) { b.setRange(2,0); }
134  BoxArray ba2d(std::move(bl2d));
135  const DistributionMapping& dm = mf.DistributionMap();
136  const int ncomp = 1;
137  IntVect ng = mf.nGrowVect(); ng[2]=0;
138 
139  m_fields[lev][1] = vars_old[Vars::yvel];
140  m_averages[lev][1] = std::make_unique<MultiFab>(ba2d,dm,ncomp,ng);
141  m_averages[lev][1]->setVal(Real(1e34));
142  if (m_rotate) {
143  m_rot_fields[lev][1] = std::make_unique<MultiFab>(ba,dm,ncomp,ng);
144  MultiFab::Copy(*m_rot_fields[lev][1],mf,0,0,1,ng);
145  } else {
146  m_rot_fields[lev][1] = nullptr;
147  }
148  }
149  { // CC vars
150  auto& mf = *Theta_prim;
151  // Create a 2D ba, dm, & ghost cells
152  const BoxArray& ba = mf.boxArray();
153  BoxList bl2d = ba.boxList();
154  for (auto& b : bl2d) { b.setRange(2,0); }
155  BoxArray ba2d(std::move(bl2d));
156  const DistributionMapping& dm = mf.DistributionMap();
157  const int ncomp = 1;
158  const int incomp = 1;
159  IntVect ng = mf.nGrowVect(); ng[2]=0;
160 
161  // Get field pointers
162  m_fields[lev][2] = Theta_prim.get();
163  m_fields[lev][3] = Qv_prim.get();
164  m_fields[lev][4] = Qr_prim.get();
165 
166  // Initialize remaining multifabs
167  for (int iavg(2); iavg < m_navg; ++iavg) {
168  m_averages[lev][iavg] = std::make_unique<MultiFab>(ba2d,dm,ncomp,ng);
169  m_averages[lev][iavg]->setVal(Real(1e34));
170  }
171 
172  // Default to dry
173  m_averages[lev][3]->setVal(0.0);
174 
175  if (m_rotate) {
176  m_rot_fields[lev][2] = std::make_unique<MultiFab>(ba,dm,ncomp,ng);
177  m_rot_fields[lev][3] = std::make_unique<MultiFab>(ba,dm,ncomp,ng);
178  MultiFab::Copy(*m_rot_fields[lev][2],*Theta_prim,0,0,1,ng);
179  if (Qv_prim) MultiFab::Copy(*m_rot_fields[lev][3],*Qv_prim,0,0,1,ng);
180  } else {
181  m_rot_fields[lev][2] = nullptr;
182  m_rot_fields[lev][3] = nullptr;
183  }
184 
185  // Default zref to 10 and fill will true values later
186  m_zref[lev] = std::make_unique<MultiFab>(ba2d,dm,1,ng);
187  m_zref[lev]->setVal(zref_default);
188 
189  if (use_terrain_fitted_coords && m_norm_vec && m_interp) {
190  m_x_pos[lev] = std::make_unique<MultiFab>(ba2d,dm,ncomp,ng);
191  m_y_pos[lev] = std::make_unique<MultiFab>(ba2d,dm,ncomp,ng);
192  m_z_pos[lev] = std::make_unique<MultiFab>(ba2d,dm,ncomp,ng);
193  } else if (use_terrain_fitted_coords && m_interp) {
194  m_x_pos[lev] = std::make_unique<MultiFab>(ba2d,dm,ncomp,ng);
195  m_y_pos[lev] = std::make_unique<MultiFab>(ba2d,dm,ncomp,ng);
196  m_z_pos[lev] = std::make_unique<MultiFab>(ba2d,dm,ncomp,ng);
197  } else if (use_terrain_fitted_coords && m_norm_vec) {
198  m_i_indx[lev] = std::make_unique<iMultiFab>(ba2d,dm,incomp,ng);
199  m_j_indx[lev] = std::make_unique<iMultiFab>(ba2d,dm,incomp,ng);
200  m_k_indx[lev] = std::make_unique<iMultiFab>(ba2d,dm,incomp,ng);
201  } else {
202  if (!use_eb) {
203  m_k_indx[lev] = std::make_unique<iMultiFab>(ba2d,dm,incomp,ng);
204  }
205  }
206  }
207  // Nodal in z (only used with terrain stress rotations)
208  m_fields[lev][5] = vars_old[Vars::zvel];
209 
210  // Setup auxiliary data for spatial configuration & policy
211  //--------------------------------------------------------
212  if (use_terrain_fitted_coords && m_norm_vec && m_interp) { // Terrain w/ norm & w/ interpolation
213  set_norm_positions_T(lev);
214  } else if (use_terrain_fitted_coords && m_interp) { // Terrain w/ interpolation
215  set_z_positions_T(lev);
216  } else if (use_terrain_fitted_coords && m_norm_vec) { // Terrain w/ norm & w/o interpolation
217  set_norm_indices_T(lev);
218  } else if (use_terrain_fitted_coords) { // Terrain
219  set_k_indices_T(lev);
220  } else if (use_eb) { // EB
221  set_z_positions_EB(lev);
222  } else { // No Terrain
223  set_k_indices_N(lev);
224  }
225 
226  // Setup normalization data for the chosen policy
227  //--------------------------------------------------------
228  switch(m_policy) {
229  case 0: // Plane average
230  set_plane_normalization(lev);
231  break;
232  case 1: // Local region/point
233  set_region_normalization(lev);
234  break;
235  case 2: // EB
236  set_eb_normalization(lev);
237  break;
238  default:
239  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(false, "Unknown policy for MOSTAverage!");
240  }
241 
242  // Set up the exponential time filtering
243  //--------------------------------------------------------
244  if (m_t_avg) {
245  // Exponential filter function
246  m_fact_old = std::exp(-one / m_time_window);
247 
248  // Enforce discrete normalization: (mfn*val_new + mfo*val_old)
249  m_fact_new = one - m_fact_old;
250 
251  // None of the averages are initialized
252  m_t_init.resize(m_maxlev,0);
253  }
254 
255  // Corrections to the mean surface velocity
256  m_Vsg = Vector<Real>(m_maxlev, zero);
257  if (include_subgrid_vel) {
258  if (include_subgrid_vel) {
259  Print() << "Subgrid velocity scale correction at level : " << lev << ' ';
260  const auto dxArr = m_geom[lev].CellSizeArray();
261  Real dx = std::sqrt(dxArr[0]*dxArr[1]);
262  if (dx > Real(5000.)) {
263  m_Vsg[lev] = Real(0.32) * std::pow(dx/Real(5000.)-1, Real(0.33));
264  }
265  Print() << m_Vsg[lev] << std::endl;
266  }
267  }
268 }
MOSTAverage::set_z_positions_EB
void set_z_positions_EB(const int &lev)
Definition: ERF_MOSTAverage.cpp:558
MOSTAverage::set_region_normalization
void set_region_normalization(const int &)
Definition: ERF_MOSTAverage.H:66
MOSTAverage::set_z_positions_T
void set_z_positions_T(const int &lev)
Definition: ERF_MOSTAverage.cpp:748
MOSTAverage::set_norm_positions_T
void set_norm_positions_T(const int &lev)
Definition: ERF_MOSTAverage.cpp:811
MOSTAverage::set_k_indices_T
void set_k_indices_T(const int &lev)
Definition: ERF_MOSTAverage.cpp:580
MOSTAverage::set_plane_normalization
void set_plane_normalization(const int &lev)
Definition: ERF_MOSTAverage.cpp:357
MOSTAverage::m_terrain_type
TerrainType m_terrain_type
Definition: ERF_MOSTAverage.H:208
MOSTAverage::m_norm_vec
bool m_norm_vec
Definition: ERF_MOSTAverage.H:241
MOSTAverage::m_nvar
int m_nvar
Definition: ERF_MOSTAverage.H:212
MOSTAverage::set_norm_indices_T
void set_norm_indices_T(const int &lev)
Definition: ERF_MOSTAverage.cpp:660
MOSTAverage::m_maxlev
int m_maxlev
Definition: ERF_MOSTAverage.H:214
MOSTAverage::set_eb_normalization
void set_eb_normalization(const int &lev)
Definition: ERF_MOSTAverage.cpp:414
MOSTAverage::include_subgrid_vel
bool include_subgrid_vel
Definition: ERF_MOSTAverage.H:257
MOSTAverage::m_time_window
amrex::Real m_time_window
Definition: ERF_MOSTAverage.H:252
MOSTAverage::set_k_indices_N
void set_k_indices_N(const int &lev)
Definition: ERF_MOSTAverage.cpp:502
MOSTAverage::zref_default
const amrex::Real zref_default
Definition: ERF_MOSTAverage.H:262
Vars::xvel
@ xvel
Definition: ERF_IndexDefines.H:175
Vars::zvel
@ zvel
Definition: ERF_IndexDefines.H:177
Vars::yvel
@ yvel
Definition: ERF_IndexDefines.H:176

Referenced by SurfaceLayer::make_SurfaceLayer_at_level().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ operator=() [1/2]

MOSTAverage& MOSTAverage::operator= ( const MOSTAverage &  other)
delete

◆ operator=() [2/2]

MOSTAverage& MOSTAverage::operator= ( MOSTAverage &&  other)
deletenoexcept

◆ set_eb_normalization()

void MOSTAverage::set_eb_normalization ( const int &  lev)

Function to compute normalization for average over EB.

415 {
416  AMREX_ALWAYS_ASSERT(m_eb_vec[lev] != nullptr);
417 
418  // Get EB data - need both cell-centered and face-centered flags
419  const auto& cc_flags = m_eb_vec[lev]->get_const_factory()->getMultiEBCellFlagFab();
420 
421  // Get area fractions for different centerings
422  auto cc_afrac = m_eb_vec[lev]->get_const_factory()->getAreaFrac(); // Cell-centered area fractions (Array of 3 MultiCutFab*)
423 
424  // Initialize storage
425  m_total_bndry_area.resize(m_maxlev);
426  m_total_bndry_area[lev].resize(m_navg, zero);
427  m_plane_average.resize(m_maxlev);
428  m_plane_average[lev].resize(m_navg, zero);
429 
430  // Compute total area for each field type based on its centering
431  // iavg: 0=U(xface), 1=V(yface), 2=T(cc), 3=Qv(cc), 4=Tv(cc), 5=Umag(cc)
432 
433  // Get geometry for cell sizes
434  auto const& dx_arr = m_geom[lev].CellSizeArray();
435  Real dx = dx_arr[0];
436  Real dy = dx_arr[1];
437  Real dz = dx_arr[2];
438 
439  // GPU array to accumulate areas
440  Gpu::DeviceVector<Real> area_vec(m_navg, zero);
441  Real* area_device = area_vec.data();
442 
443  // All fields are now averaged on cell-centered grid
444  // Compute total area once on cell-centered grid and use for all iavg
445  Real total_area = zero;
446 #ifdef _OPENMP
447 #pragma omp parallel if (Gpu::notInLaunchRegion())
448 #endif
449  for (MFIter mfi(cc_flags, TileNoZ()); mfi.isValid(); ++mfi) {
450  const auto& flag = cc_flags[mfi];
451 
452  // Skip boxes that are not singlevalued (MultiCutFab only has data for singlevalued boxes)
453  if (flag.getType() != FabType::singlevalued) continue;
454 
455  Box bx = mfi.tilebox();
456  auto const flag_arr = flag.const_array();
457  auto const afrac_x = cc_afrac[0]->const_array(mfi);
458  auto const afrac_y = cc_afrac[1]->const_array(mfi);
459  auto const afrac_z = cc_afrac[2]->const_array(mfi);
460 
461  // Sum area only for cut cells using atomic reduction
462  ParallelFor(Gpu::KernelInfo().setReduction(true), bx, [=]
463  AMREX_GPU_DEVICE(int i, int j, int k, Gpu::Handler const& handler) noexcept
464  {
465  if (flag_arr(i,j,k).isSingleValued()) {
466  // Compute area from face-centered area fractions
467  Real axm = afrac_x(i ,j ,k );
468  Real axp = afrac_x(i+1,j ,k );
469  Real aym = afrac_y(i ,j ,k );
470  Real ayp = afrac_y(i ,j+1,k );
471  Real azm = afrac_z(i ,j ,k );
472  Real azp = afrac_z(i ,j ,k+1);
473 
474  Real adx = (axm - axp) * dy * dz;
475  Real ady = (aym - ayp) * dx * dz;
476  Real adz = (azm - azp) * dx * dy;
477 
478  Real area = std::sqrt(adx*adx + ady*ady + adz*adz);
479 
480  Gpu::deviceReduceSum(&area_device[0], area, handler);
481  }
482  });
483  }
484 
485  // Copy to host and sum across MPI ranks
486  Gpu::copy(Gpu::deviceToHost, area_vec.begin(), area_vec.begin() + 1, &total_area);
487  ParallelDescriptor::ReduceRealSum(&total_area, 1);
488 
489  // Set the same total area for all iavg
490  for (int iavg = 0; iavg < m_navg; ++iavg) {
491  m_total_bndry_area[lev][iavg] = total_area;
492  }
493 
494  Print() << "EB surface area on cell-centerd grid at level " << lev << ": " << total_area << std::endl;
495 }

Referenced by make_MOSTAverage_at_level().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ set_k_indices_N()

void MOSTAverage::set_k_indices_N ( const int &  lev)

Function to set K indices without terrain.

503 {
504  ParmParse pp(m_pp_prefix);
505  Real zref_tmp = zref_default;
506  auto read_z = pp.query("most.zref",zref_tmp);
507  auto read_k = pp.queryarr("most.k_arr_in",m_k_in);
508 
509  // Default behavior is to use the first cell center
510  if (!read_z && !read_k) {
511  Real m_zlo = m_geom[0].ProbLo(2);
512  Real m_dz = m_geom[0].CellSize(2);
513  zref_tmp = m_zlo + myhalf * m_dz;
514  m_zref[lev]->setVal( zref_tmp );
515  Print() << "Reference height for MOST set to " << zref_tmp << std::endl;
516  read_z = true;
517  }
518 
519  // Specify z_ref & compute k_indx (z_ref takes precedence)
520  if (read_z) {
521  Real m_zlo = m_geom[lev].ProbLo(2);
522  Real m_zhi = m_geom[lev].ProbHi(2);
523  Real m_dz = m_geom[lev].CellSize(2);
524 
525  amrex::ignore_unused(m_zhi);
526 
527  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(zref_tmp >= m_zlo + myhalf * m_dz,
528  "Query point must be past first z-cell!");
529 
530  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(zref_tmp <= m_zhi - myhalf * m_dz,
531  "Query point must be below the last z-cell!");
532 
533  int lk = static_cast<int>(floor((zref_tmp - m_zlo) / m_dz - myhalf));
534 
535  m_zref[lev]->setVal( (lk + myhalf) * m_dz + m_zlo );
536 
537  AMREX_ALWAYS_ASSERT(lk >= m_radius);
538 
539  m_k_indx[lev]->setVal(lk);
540  // Specified k_indx & compute z_ref
541  } else if (read_k) {
542  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_k_in[lev] >= m_radius,
543  "K index must be larger than averaging radius!");
544  m_k_indx[lev]->setVal(m_k_in[lev]);
545 
546  // TODO: check that z_ref is constant across levels
547  Real m_zlo = m_geom[0].ProbLo(2);
548  Real m_dz = m_geom[0].CellSize(2);
549  m_zref[lev]->setVal( ((Real)m_k_in[0] + myhalf) * m_dz + m_zlo );
550  }
551 }
pp
ParmParse pp("prob")
MOSTAverage::m_pp_prefix
std::string m_pp_prefix
Definition: ERF_MOSTAverage.H:206
MOSTAverage::m_k_in
amrex::Vector< int > m_k_in
Definition: ERF_MOSTAverage.H:236

Referenced by make_MOSTAverage_at_level().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ set_k_indices_T()

void MOSTAverage::set_k_indices_T ( const int &  lev)

Function to set K indices with terrain (w/o terrain normals or interpolation).

581 {
582  // Peel back the level
583  auto& fields = m_fields[lev];
584 
585  // MFIter over CC data
586  int imf_cc = 2;
587 
588  ParmParse pp(m_pp_prefix);
589  Real zref_tmp = zref_default;
590  auto read_z = pp.query("most.zref",zref_tmp);
591  auto read_k = pp.queryarr("most.k_arr_in",m_k_in);
592  int klo = m_geom[lev].Domain().smallEnd(2);
593 
594  // Allow default zref
595  if (!read_z) {
596  Print() << "most.zref not specified, query distance default is " << zref_tmp << std::endl;
597  read_z = true;
598  }
599 
600  // No default behavior with terrain (we can't tell the difference between
601  // vertical grid stretching and true terrain)
602  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(read_z != read_k,
603  "Need to specify zref or k_arr_in for MOST");
604 
605  // Capture for device
606  Real d_zref = zref_tmp;
607  Real d_radius = static_cast<Real>(m_radius);
608  amrex::ignore_unused(d_radius);
609 
610  // Specify z_ref & compute k_indx (z_ref takes precedence)
611  if (read_z) {
612  int kmax = m_geom[lev].Domain().bigEnd(2);
613  for (MFIter mfi(*fields[imf_cc], TileNoZ()); mfi.isValid(); ++mfi) {
614  Box npbx = mfi.tilebox(IntVect(1,1,0),IntVect(1,1,0));
615 
616  if (npbx.smallEnd(2) != klo) { continue; }
617 
618  npbx.makeSlab(2,klo);
619 
620  const auto z_phys_arr = m_z_phys_nd[lev]->const_array(mfi);
621  auto k_arr = m_k_indx[lev]->array(mfi);
622  auto zref_arr = m_zref[lev]->array(mfi);
623  ParallelFor(npbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
624  {
625  k_arr(i,j,k) = klo;
626  bool found = false;
627  Real z_bot_face = fourth * ( z_phys_arr(i ,j ,k) + z_phys_arr(i+1,j ,k)
628  + z_phys_arr(i ,j+1,k) + z_phys_arr(i+1,j+1,k) );
629  Real z_target = z_bot_face + d_zref;
630  for (int lk(klo); lk<=kmax; ++lk) {
631  Real z_lo = fourth * ( z_phys_arr(i,j ,lk ) + z_phys_arr(i+1,j ,lk )
632  + z_phys_arr(i,j+1,lk ) + z_phys_arr(i+1,j+1,lk ) );
633  Real z_hi = fourth * ( z_phys_arr(i,j ,lk+1) + z_phys_arr(i+1,j ,lk+1)
634  + z_phys_arr(i,j+1,lk+1) + z_phys_arr(i+1,j+1,lk+1) );
635  if (z_target > z_lo && z_target < z_hi){
636  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(lk >= d_radius,
637  "K index must be larger than averaging radius!");
638  k_arr(i,j,0) = lk;
639  zref_arr(i,j,0) = myhalf * (z_hi + z_lo) - z_bot_face;
640  found = true;
641  break;
642  }
643  }
644  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(found,
645  "zref not found with terrain!");
646  });
647  }
648  // Specified k_indx & compute z_ref
649  } else if (read_k) {
650  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(false, "Specified k-indx with terrain not implemented!");
651  }
652 }
fourth
constexpr amrex::Real fourth
Definition: ERF_Constants.H:12

Referenced by make_MOSTAverage_at_level().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ set_norm_indices_T()

void MOSTAverage::set_norm_indices_T ( const int &  lev)

Function to set I,J,K indices with terrain normals (w/o interpolation).

661 {
662  // Peel back the level
663  auto& fields = m_fields[lev];
664 
665  // MFIter over CC data
666  int imf_cc = 2;
667 
668  ParmParse pp(m_pp_prefix);
669  Real zref_tmp = zref_default;
670  auto read_zref = pp.query("most.zref",zref_tmp);
671  if (!read_zref) {
672  Print() << "most.zref not specified, query distance default is " << zref_tmp << std::endl;
673  }
674  int klo = m_geom[lev].Domain().smallEnd(2);
675 
676  // Capture for device
677  Real d_zref = zref_tmp;
678  Real d_radius = static_cast<Real>(m_radius);
679 
680  const auto dxInv = m_geom[lev].InvCellSizeArray();
681  IntVect ng = m_k_indx[lev]->nGrowVect(); ng[2]=0;
682  for (MFIter mfi(*fields[imf_cc], TileNoZ()); mfi.isValid(); ++mfi) {
683  Box npbx = mfi.tilebox(IntVect(1,1,0),IntVect(1,1,0));
684 
685  if (npbx.smallEnd(2) != klo) { continue; }
686 
687  int kmax = npbx.bigEnd(2);
688 
689  npbx.makeSlab(2,klo);
690 
691  const auto z_phys_arr = m_z_phys_nd[lev]->const_array(mfi);
692  auto i_arr = m_i_indx[lev]->array(mfi);
693  auto j_arr = m_j_indx[lev]->array(mfi);
694  auto k_arr = m_k_indx[lev]->array(mfi);
695  auto zref_arr = m_zref[lev]->array(mfi);
696  ParallelFor(npbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
697  {
698  // Elements of normal vector
699  Real met_h_xi = Compute_h_xi_AtCellCenter (i,j,k,dxInv,z_phys_arr);
700  Real met_h_eta = Compute_h_eta_AtCellCenter(i,j,k,dxInv,z_phys_arr);
701  Real mag = std::sqrt(met_h_xi*met_h_xi + met_h_eta*met_h_eta + one);
702 
703  // Unit-normal vector scaled by z_ref
704  Real delta_x = -met_h_xi/mag * d_zref;
705  Real delta_y = -met_h_eta/mag * d_zref;
706  Real delta_z = one/mag * d_zref;
707 
708  // Compute i & j as displacements (no grid stretching)
709  int delta_i = static_cast<int>(std::round(delta_x*dxInv[0]));
710  int delta_j = static_cast<int>(std::round(delta_y*dxInv[1]));
711  int i_new = i + delta_i;
712  int j_new = j + delta_j;
713  i_arr(i,j,0) = i_new;
714  j_arr(i,j,0) = j_new;
715 
716  // Search for k (grid is stretched in z)
717  Real z_bot_face = fourth * ( z_phys_arr(i ,j ,k) + z_phys_arr(i+1,j ,k)
718  + z_phys_arr(i ,j+1,k) + z_phys_arr(i+1,j+1,k) );
719  Real z_target = z_bot_face + delta_z;
720  k_arr(i,j,0) = klo;
721  zref_arr(i,j,0) = myhalf * z_bot_face +
722  Real(0.125) * ( z_phys_arr(i ,j ,k+1) + z_phys_arr(i+1,j ,k+1)
723  + z_phys_arr(i ,j+1,k+1) + z_phys_arr(i+1,j+1,k+1) );
724  for (int lk(klo); lk<=kmax; ++lk) {
725  Real z_lo = fourth * ( z_phys_arr(i_new,j_new ,lk ) + z_phys_arr(i_new+1,j_new ,lk )
726  + z_phys_arr(i_new,j_new+1,lk ) + z_phys_arr(i_new+1,j_new+1,lk ) );
727  Real z_hi = fourth * ( z_phys_arr(i_new,j_new ,lk+1) + z_phys_arr(i_new+1,j_new ,lk+1)
728  + z_phys_arr(i_new,j_new+1,lk+1) + z_phys_arr(i_new+1,j_new+1,lk+1) );
729  if (z_target > z_lo && z_target < z_hi){
730  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(lk >= d_radius,
731  "K index must be larger than averaging radius!");
732  amrex::ignore_unused(d_radius);
733  k_arr(i,j,0) = lk;
734  zref_arr(i,j,0) = myhalf * (z_hi + z_lo) - z_bot_face;
735  break;
736  }
737  }
738  });
739  }
740 }
Compute_h_eta_AtCellCenter
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_eta_AtCellCenter(const int &i, const int &j, const int &k, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &z_nd)
Definition: ERF_TerrainMetrics.H:85
Compute_h_xi_AtCellCenter
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_xi_AtCellCenter(const int &i, const int &j, const int &k, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &z_nd)
Definition: ERF_TerrainMetrics.H:70

Referenced by make_MOSTAverage_at_level().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ set_norm_positions_T()

void MOSTAverage::set_norm_positions_T ( const int &  lev)

Function to set positions with terrain and normal vector (with interpolation).

812 {
813  // Peel back the level
814  auto& fields = m_fields[lev];
815 
816  // MFIter over CC data
817  int imf_cc = 2;
818 
819  ParmParse pp(m_pp_prefix);
820  Real zref_tmp = zref_default;
821  auto read_zref = pp.query("most.zref",zref_tmp);
822  if (!read_zref) {
823  Print() << "most.zref not specified, query distance default is " << zref_tmp << std::endl;
824  }
825  int klo = m_geom[lev].Domain().smallEnd(2);
826 
827  // Capture for device
828  Real d_zref = zref_tmp;
829  const auto plo = m_geom[lev].ProbLoArray();
830 
831  RealVect base;
832  const auto dx = m_geom[lev].CellSizeArray();
833  const auto dxInv = m_geom[lev].InvCellSizeArray();
834  IntVect ng = m_x_pos[lev]->nGrowVect(); ng[2]=0;
835  for (MFIter mfi(*fields[imf_cc], TileNoZ()); mfi.isValid(); ++mfi) {
836  Box npbx = mfi.tilebox(IntVect(1,1,0),IntVect(1,1,0));
837  Box gtbx = mfi.growntilebox(ng);
838  RealBox grb{gtbx,dx.data(),base.dataPtr()};
839 
840  if (npbx.smallEnd(2) != klo) { continue; }
841 
842  npbx.makeSlab(2,klo);
843 
844  const auto z_phys_arr = m_z_phys_nd[lev]->const_array(mfi);
845  auto x_pos_arr = m_x_pos[lev]->array(mfi);
846  auto y_pos_arr = m_y_pos[lev]->array(mfi);
847  auto z_pos_arr = m_z_pos[lev]->array(mfi);
848  auto zref_arr = m_zref[lev]->array(mfi);
849  ParallelFor(npbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
850  {
851  // Elements of normal vector
852  Real met_h_xi = Compute_h_xi_AtCellCenter (i,j,k,dxInv,z_phys_arr);
853  Real met_h_eta = Compute_h_eta_AtCellCenter(i,j,k,dxInv,z_phys_arr);
854  Real imag = one / std::sqrt(met_h_xi*met_h_xi + met_h_eta*met_h_eta + one);
855 
856  // Unit-normal vector scaled by z_ref
857  Real delta_x = -met_h_xi * imag * d_zref;
858  Real delta_y = -met_h_eta * imag * d_zref;
859  Real delta_z = imag * d_zref;
860 
861  // Position of the current node (indx:0,0,1)
862  Real x0 = plo[0] + ((Real) i + myhalf) * dx[0];
863  Real y0 = plo[1] + ((Real) j + myhalf) * dx[1];
864 
865  // Final position at end of vector
866  x_pos_arr(i,j,0) = x0 + delta_x;
867  y_pos_arr(i,j,0) = y0 + delta_y;
868  Real z_bot_face = fourth * ( z_phys_arr(i ,j ,k) + z_phys_arr(i+1,j ,k)
869  + z_phys_arr(i ,j+1,k) + z_phys_arr(i+1,j+1,k) );
870  z_pos_arr(i,j,0) = z_bot_face + delta_z;
871 
872  // NOTE: Normal vector end point can be below the surface for concave regions.
873  // Here we protect against that by augmenting the normal if needed.
874  int i_new = (int) ((x_pos_arr(i,j,0) - plo[0]) / dx[0] - myhalf);
875  int j_new = (int) ((y_pos_arr(i,j,0) - plo[1]) / dx[1] - myhalf);
876  Real z_new_bot_face = fourth * ( z_phys_arr(i_new,j_new ,k) + z_phys_arr(i_new+1,j_new ,k)
877  + z_phys_arr(i_new,j_new+1,k) + z_phys_arr(i_new+1,j_new+1,k) );
878  if (z_pos_arr(i,j,0) < z_new_bot_face) {
879  z_pos_arr(i,j,0) = z_new_bot_face + delta_z;
880  }
881 
882  zref_arr(i,j,0) = delta_z;
883 
884  // Destination position must be contained on the current process!
885  Real pos[] = {x_pos_arr(i,j,0)-plo[0],y_pos_arr(i,j,0)-plo[1],myhalf*dx[2]};
886  amrex::ignore_unused(pos);
887  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(grb.contains(&pos[0]),
888  "Query point outside of proc domain!");
889  });
890  }
891 }

Referenced by make_MOSTAverage_at_level().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ set_plane_normalization()

void MOSTAverage::set_plane_normalization ( const int &  lev)

Function to compute normalization for plane average.

358 {
359  // Cells per plane and temp avg storage
360  m_ncell_plane.resize(m_maxlev);
361  m_plane_average.resize(m_maxlev);
362 
363  // True domain not used for normalization
364  Box domain = m_geom[lev].Domain();
365 
366  // NOTE: Level 0 spans the whole domain, but finer
367  // levels do not have such a restriction.
368  // For now, use the bounding box of the boxArray
369  // for normalization, consistent with avg routine.
370 
371  // Bounded box of CC data used for normalization
372  Box bnd_bx = (m_fields[lev][2]->boxArray()).minimalBox();
373 
374  // NOTE: Bounding box must lie on the periodic boundaries
375  // in order to trip the is_per flag
376 
377  // Num components, plane avg, cells per plane
378  Array<int,AMREX_SPACEDIM> is_per = {0,0,0};
379  for (int idim(0); idim < AMREX_SPACEDIM-1; ++idim) {
380  if ( m_geom[lev].isPeriodic(idim) &&
381  bnd_bx.bigEnd(idim)==domain.bigEnd(idim) &&
382  bnd_bx.smallEnd(idim)==domain.smallEnd(idim) ) { is_per[idim] = 1; }
383  }
384 
385  m_ncell_plane[lev].resize(m_navg);
386  m_plane_average[lev].resize(m_navg);
387  for (int iavg(0); iavg < m_navg; ++iavg) {
388  // Convert bnd_bx to current index type
389  IndexType ixt = m_averages[lev][iavg]->boxArray().ixType();
390  bnd_bx.convert(ixt);
391  IntVect bnd_bx_lo(bnd_bx.loVect());
392  IntVect bnd_bx_hi(bnd_bx.hiVect());
393 
394  m_plane_average[lev][iavg] = zero;
395 
396  m_ncell_plane[lev][iavg] = 1;
397  for (int idim(0); idim < AMREX_SPACEDIM; ++idim) {
398  if (idim != 2) {
399  if (ixt.nodeCentered(idim) && is_per[idim]) {
400  m_ncell_plane[lev][iavg] *= (bnd_bx_hi[idim] - bnd_bx_lo[idim]);
401  } else {
402  m_ncell_plane[lev][iavg] *= (bnd_bx_hi[idim] - bnd_bx_lo[idim] + 1);
403  }
404  }
405  } // idim
406  } // iavg
407 }

Referenced by make_MOSTAverage_at_level().

Here is the caller graph for this function:

◆ set_region_normalization()

void MOSTAverage::set_region_normalization ( const int &  )
inline
67  {m_ncell_region = (2 * m_radius + 1) * (2 * m_radius + 1) * (2 * m_radius + 1);}

Referenced by make_MOSTAverage_at_level().

Here is the caller graph for this function:

◆ set_rotated_fields()

void MOSTAverage::set_rotated_fields ( const int &  lev)

Function to set the rotated velocities.

299 {
300  // Peel back the level
301  auto& fields = m_fields[lev];
302  auto& rot_fields = m_rot_fields[lev];
303  auto z_phys_nd = m_z_phys_nd[lev];
304 
305  // Inverse grid size
306  const auto dxInv = m_geom[lev].InvCellSizeArray();
307 
308  // Single MFIter over CC data
309  int imf_cc = 2;
310 
311  // Populate rotated U & V for terrain
312 #ifdef _OPENMP
313 #pragma omp parallel if (Gpu::notInLaunchRegion())
314 #endif
315  for (MFIter mfi(*fields[imf_cc], TileNoZ()); mfi.isValid(); ++mfi) {
316  Box ubx = mfi.tilebox(IntVect(1,0,0));
317  Box vbx = mfi.tilebox(IntVect(0,1,0));
318 
319  const Array4<const Real>& z_phys_arr = z_phys_nd->const_array(mfi);
320 
321  const Array4<const Real>& u_arr = fields[0]->const_array(mfi);
322  const Array4<const Real>& v_arr = fields[1]->const_array(mfi);
323  const Array4<const Real>& w_arr = fields[5]->const_array(mfi);
324 
325  const Array4<Real>& u_rot_arr = rot_fields[0]->array(mfi);
326  const Array4<Real>& v_rot_arr = rot_fields[1]->array(mfi);
327 
328  // U rotated magnitude
329  ParallelFor(ubx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
330  {
331  // Elements of first tangent vector
332  Real met_h_xi = Compute_h_xi_AtIface(i,j,k,dxInv,z_phys_arr);
333  u_rot_arr(i,j,k) = (u_arr(i,j,k) + met_h_xi*w_arr(i,j,k))
334  / std::sqrt(met_h_xi*met_h_xi + one);
335  });
336 
337  // V rotated magnitude
338  ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
339  {
340  // Elements of second tangent vector
341  Real met_h_eta = Compute_h_eta_AtJface(i,j,k,dxInv,z_phys_arr);
342  v_rot_arr(i,j,k) = (v_arr(i,j,k) + met_h_eta*w_arr(i,j,k))
343  / std::sqrt(met_h_eta*met_h_eta + one);
344  });
345  }
346 
347  // Direct copy of other scalar variables
348  MultiFab::Copy(*rot_fields[2],*fields[2],0,0,1,rot_fields[2]->nGrowVect());
349  if (fields[3]) MultiFab::Copy(*rot_fields[3],*fields[3],0,0,1,rot_fields[3]->nGrowVect());
350 }
Compute_h_xi_AtIface
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_xi_AtIface(const int &i, const int &j, const int &k, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &z_nd)
Definition: ERF_TerrainMetrics.H:117
Compute_h_eta_AtJface
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_eta_AtJface(const int &i, const int &j, const int &k, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &z_nd)
Definition: ERF_TerrainMetrics.H:170

Referenced by compute_averages().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ set_z_positions_EB()

void MOSTAverage::set_z_positions_EB ( const int &  lev)

Function to set K indices for EB.

559 {
560  Real zref_tmp = zref_default;
561  ParmParse pp(m_pp_prefix);
562  auto read_z = pp.query("most.zref",zref_tmp);
563 
564  if (!read_z) {
565  m_zref[lev]->setVal( zref_tmp );
566  // Default behavior is to use the first cell center
567  } else {
568  Real m_dz = m_geom[0].CellSize(2);
569  zref_tmp = myhalf * m_dz;
570  m_zref[lev]->setVal( zref_tmp );
571  Print() << "Reference height for MOST set to " << zref_tmp << std::endl;
572  }
573 }

Referenced by make_MOSTAverage_at_level().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ set_z_positions_T()

void MOSTAverage::set_z_positions_T ( const int &  lev)

Function to set positions with terrain and e_z vector (with interpolation but no terrain normals)

749 {
750  // Peel back the level
751  auto& fields = m_fields[lev];
752 
753  // MFIter over CC data
754  int imf_cc = 2;
755 
756  ParmParse pp(m_pp_prefix);
757  Real zref_tmp = zref_default;
758  auto read_zref = pp.query("most.zref",zref_tmp);
759  if (!read_zref) {
760  Print() << "most.zref not specified, query distance default is " << zref_tmp << std::endl;
761  } else {
762  m_zref[lev]->setVal(zref_tmp);
763  }
764  int klo = m_geom[lev].Domain().smallEnd(2);
765 
766  // Capture for device
767  Real d_zref = zref_tmp;
768  const auto plo = m_geom[lev].ProbLoArray();
769 
770  RealVect base;
771  const auto dx = m_geom[lev].CellSizeArray();
772  IntVect ng = m_x_pos[lev]->nGrowVect(); ng[2]=0;
773  for (MFIter mfi(*fields[imf_cc], TileNoZ()); mfi.isValid(); ++mfi) {
774  Box npbx = mfi.tilebox(IntVect(1,1,0),IntVect(1,1,0));
775  Box gtbx = mfi.growntilebox(ng);
776 
777  if (npbx.smallEnd(2) != klo) { continue; }
778 
779  npbx.makeSlab(2,klo);
780 
781  RealBox grb{gtbx,dx.data(),base.dataPtr()};
782 
783  const auto z_phys_arr = m_z_phys_nd[lev]->const_array(mfi);
784  auto x_pos_arr = m_x_pos[lev]->array(mfi);
785  auto y_pos_arr = m_y_pos[lev]->array(mfi);
786  auto z_pos_arr = m_z_pos[lev]->array(mfi);
787  ParallelFor(npbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
788  {
789  // Final position at end of vector
790  x_pos_arr(i,j,0) = plo[0] + ((Real) i + myhalf) * dx[0];
791  y_pos_arr(i,j,0) = plo[1] + ((Real) j + myhalf) * dx[1];
792  Real z_bot_face = fourth * ( z_phys_arr(i ,j ,k) + z_phys_arr(i+1,j ,k)
793  + z_phys_arr(i ,j+1,k) + z_phys_arr(i+1,j+1,k) );
794  z_pos_arr(i,j,0) = z_bot_face + d_zref;
795 
796  // Destination position must be contained on the current process!
797  Real pos[] = {x_pos_arr(i,j,0)-plo[0],y_pos_arr(i,j,0)-plo[1],myhalf*dx[2]};
798  amrex::ignore_unused(pos);
799  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(grb.contains(&pos[0]),
800  "Query point outside of proc domain!");
801  });
802  }
803 }

Referenced by make_MOSTAverage_at_level().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ trilinear_interp_T()

AMREX_GPU_HOST_DEVICE static AMREX_INLINE void MOSTAverage::trilinear_interp_T ( const amrex::Real &  xp,
const amrex::Real &  yp,
const amrex::Real &  zp,
amrex::Real *  interp_vals,
amrex::Array4< amrex::Real const > const &  interp_array,
amrex::Array4< amrex::Real const > const &  z_arr,
const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &  plo,
const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &  dxi,
const int  interp_comp 
)
inlinestatic

Function to compute trilinear interpolation with terrain.

Parameters
[in]xpX-position
[in]ypY-position
[out]interp_valsValues interpolated
[in]interp_arrayArray to interpolate on
[in]z_arrPhysical heights
[in]ploProblem lower bounds
[in]dxiInverse cell size array
[in]interp_compNumber of components to interpolate
143  {
144  // Search to get z/k
145  bool found = false;
146  int kmax = ubound(z_arr).z;
147  amrex::Real zval = zero;
148  amrex::Real z_target = zp;
149 
150  // Map position to i,j (must be same mapping in cpp file)
151  amrex::Real ireal = (xp - plo[0]) * dxi[0];
152  amrex::Real jreal = (yp - plo[1]) * dxi[1];
153  int i_new = (int) (ireal - myhalf);
154  int j_new = (int) (jreal - myhalf);
155 
156  for (int lk(0); lk<kmax; ++lk) {
157  amrex::Real z_lo = fourth * ( z_arr(i_new,j_new ,lk ) + z_arr(i_new+1,j_new ,lk )
158  + z_arr(i_new,j_new+1,lk ) + z_arr(i_new+1,j_new+1,lk ) );
159  amrex::Real z_hi = fourth * ( z_arr(i_new,j_new ,lk+1) + z_arr(i_new+1,j_new ,lk+1)
160  + z_arr(i_new,j_new+1,lk+1) + z_arr(i_new+1,j_new+1,lk+1) );
161  if (z_target > z_lo && z_target < z_hi){
162  found = true;
163  zval = (amrex::Real) lk + ((z_target - z_lo) / (z_hi - z_lo)) + myhalf;
164  break;
165  }
166  }
167 
168  amrex::ignore_unused(found);
169  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(found, "MOSTAverage: Height above terrain not found, try increasing z_ref!");
170 
171  // NOTE: This is the point ahead of the current i,j (e.g. i/j_new + 1)
172  const amrex::RealVect lx(ireal + myhalf, jreal + myhalf, zval);
173 
174  const amrex::IntVect ijk = lx.floor();
175 
176  int i = ijk[0]; int j = ijk[1]; int k = ijk[2];
177 
178  // Convert ijk (IntVect) to a RealVect explicitly
179  amrex::RealVect ijk_r(static_cast<amrex::Real>(ijk[0]),
180  static_cast<amrex::Real>(ijk[1]),
181  static_cast<amrex::Real>(ijk[2]));
182 
183  // Weights
184  const amrex::RealVect sx_hi = lx - ijk_r;
185  const amrex::RealVect sx_lo = one - sx_hi;
186 
187  for (int n = 0; n < interp_comp; n++) {
188  interp_vals[n] = sx_lo[0]*sx_lo[1]*sx_lo[2]*interp_array(i-1, j-1, k-1,n) +
189  sx_lo[0]*sx_lo[1]*sx_hi[2]*interp_array(i-1, j-1, k ,n) +
190  sx_lo[0]*sx_hi[1]*sx_lo[2]*interp_array(i-1, j , k-1,n) +
191  sx_lo[0]*sx_hi[1]*sx_hi[2]*interp_array(i-1, j , k ,n) +
192  sx_hi[0]*sx_lo[1]*sx_lo[2]*interp_array(i , j-1, k-1,n) +
193  sx_hi[0]*sx_lo[1]*sx_hi[2]*interp_array(i , j-1, k ,n) +
194  sx_hi[0]*sx_hi[1]*sx_lo[2]*interp_array(i , j , k-1,n) +
195  sx_hi[0]*sx_hi[1]*sx_hi[2]*interp_array(i , j , k ,n);
196  }
197  }

Referenced by compute_plane_averages(), and compute_region_averages().

Here is the caller graph for this function:

◆ update_field_ptrs()

void MOSTAverage::update_field_ptrs ( const int &  lev,
amrex::Vector< amrex::Vector< amrex::MultiFab >> &  vars_old,
amrex::Vector< std::unique_ptr< amrex::MultiFab >> &  Theta_prim,
amrex::Vector< std::unique_ptr< amrex::MultiFab >> &  Qv_prim,
amrex::Vector< std::unique_ptr< amrex::MultiFab >> &  Qr_prim 
)

Function to reset the pointers to field variables.

Parameters
[in]levCurrent level
[in]vars_oldConserved variables at each level
[in]Theta_primPrimitive theta component at each level
284 {
285  m_fields[lev][0] = &vars_old[lev][Vars::xvel];
286  m_fields[lev][1] = &vars_old[lev][Vars::yvel];
287  m_fields[lev][2] = Theta_prim[lev].get();
288  m_fields[lev][3] = Qv_prim[lev].get();
289  m_fields[lev][4] = Qr_prim[lev].get();
290  m_fields[lev][5] = &vars_old[lev][Vars::zvel];
291 }

Referenced by SurfaceLayer::update_mac_ptrs().

Here is the caller graph for this function:

◆ write_averages()

void MOSTAverage::write_averages ( const int &  lev)

Function to write averages to text file.

Parameters
[in]levCurrent level
2040 {
2041  // Peel back the level
2042  auto& fields = m_fields[lev];
2043  auto& averages = m_averages[lev];
2044 
2045  // MFIter on CC
2046  int imf_cc = 2;
2047 
2048  int klo = m_geom[lev].Domain().smallEnd(2);
2049 
2050  int navg = m_navg - 1;
2051 
2052  std::ofstream ofile;
2053  ofile.open ("MOST_averages.txt");
2054  ofile << "Averages computed via MOSTAverages class:\n";
2055 
2056  for (MFIter mfi(*fields[imf_cc], TileNoZ()); mfi.isValid(); ++mfi) {
2057  Box pbx = mfi.tilebox();
2058 
2059  if(pbx.smallEnd(2) != klo) { continue; }
2060 
2061  pbx.makeSlab(2,klo);
2062 
2063  int il = pbx.smallEnd(0); int iu = pbx.bigEnd(0);
2064  int jl = pbx.smallEnd(1); int ju = pbx.bigEnd(1);
2065 
2066  for (int j(jl); j <= ju; ++j) {
2067  for (int i(il); i <= iu; ++i) {
2068  ofile << "(I,J): " << "(" << i << "," << j << ")" << "\n";
2069  int k = 0;
2070  for (int iavg(0); iavg <= navg; ++iavg) {
2071  auto mf_arr = averages[iavg]->array(mfi);
2072  ofile << "iavg val: "
2073  << iavg << ' '
2074  << mf_arr(i,j,k) << "\n";
2075  }
2076  ofile << "\n";
2077  }
2078  }
2079  }
2080  ofile.close();
2081 }
Here is the call graph for this function:

◆ write_k_indices()

void MOSTAverage::write_k_indices ( const int &  lev)

Function to write the K indices to text file.

Parameters
[in]levCurrent level
1894 {
1895  // Peel back the level
1896  auto& fields = m_fields[lev];
1897  auto& k_indx = m_k_indx[lev];
1898 
1899  // MFIter on CC
1900  int imf_cc = 2;
1901 
1902  int klo = m_geom[lev].Domain().smallEnd(2);
1903 
1904  std::ofstream ofile;
1905  ofile.open ("MOST_k_indices.txt");
1906  ofile << "K indices used to compute averages via MOSTAverages class:\n";
1907 
1908  for (MFIter mfi(*fields[imf_cc], TileNoZ()); mfi.isValid(); ++mfi) {
1909  Box pbx = mfi.tilebox();
1910 
1911  if(pbx.smallEnd(2) != klo) { continue; }
1912 
1913  pbx.makeSlab(2,klo);
1914 
1915  int il = pbx.smallEnd(0); int iu = pbx.bigEnd(0);
1916  int jl = pbx.smallEnd(1); int ju = pbx.bigEnd(1);
1917 
1918  auto k_arr = k_indx->array(mfi);
1919 
1920  for (int j(jl); j <= ju; ++j) {
1921  for (int i(il); i <= iu; ++i) {
1922  ofile << "(I,J): " << "(" << i << "," << j << ")" << "\n";
1923  int k = 0;
1924  ofile << "K_ind: "
1925  << k_arr(i,j,k) << "\n";
1926  ofile << "\n";
1927  }
1928  }
1929  }
1930  ofile.close();
1931 }
Here is the call graph for this function:

◆ write_norm_indices()

void MOSTAverage::write_norm_indices ( const int &  lev)

Function to write I,J,K indices to text file.

Parameters
[in]levCurrent level
1941 {
1942  // Peel back the level
1943  auto& fields = m_fields[lev];
1944  auto& k_indx = m_k_indx[lev];
1945  auto& j_indx = m_j_indx[lev];
1946  auto& i_indx = m_i_indx[lev];
1947 
1948  // MFIter on CC
1949  int imf_cc = 2;
1950 
1951  int klo = m_geom[lev].Domain().smallEnd(2);
1952 
1953  std::ofstream ofile;
1954  ofile.open ("MOST_ijk_indices.txt");
1955  ofile << "IJK indices used to compute averages via MOSTAverages class:\n";
1956 
1957  for (MFIter mfi(*fields[imf_cc], TileNoZ()); mfi.isValid(); ++mfi) {
1958  Box pbx = mfi.tilebox();
1959 
1960  if(pbx.smallEnd(2) != klo) { continue; }
1961 
1962  pbx.makeSlab(2,klo);
1963 
1964  int il = pbx.smallEnd(0); int iu = pbx.bigEnd(0);
1965  int jl = pbx.smallEnd(1); int ju = pbx.bigEnd(1);
1966 
1967  auto k_arr = k_indx->array(mfi);
1968  auto j_arr = j_indx ? j_indx->array(mfi) : Array4<int> {};
1969  auto i_arr = i_indx ? i_indx->array(mfi) : Array4<int> {};
1970 
1971  for (int j(jl); j <= ju; ++j) {
1972  for (int i(il); i <= iu; ++i) {
1973  ofile << "(I1,J1,K1): " << "(" << i << "," << j << "," << 0 << ")" << "\n";
1974 
1975  int k = 0;
1976  int km = k_arr(i,j,k);
1977  int jm = j_arr ? j_arr(i,j,k) : j;
1978  int im = i_arr ? i_arr(i,j,k) : i;
1979 
1980  ofile << "(I2,J2,K2): "
1981  << "(" << im << "," << jm << "," << km << ")" << "\n";
1982  ofile << "\n";
1983  }
1984  }
1985  }
1986  ofile.close();
1987 }
Here is the call graph for this function:

◆ write_xz_positions()

void MOSTAverage::write_xz_positions ( const int &  lev,
const int &  j 
)

Function to write X & Z positions to text file.

Parameters
[in]levCurrent level
[in]jIndex in y-dir
1999 {
2000  // Peel back the level
2001  auto& fields = m_fields[lev];
2002  auto& x_pos_mf = m_x_pos[lev];
2003  auto& z_pos_mf = m_z_pos[lev];
2004 
2005  // MFIter on CC
2006  int imf_cc = 2;
2007 
2008  int klo = m_geom[lev].Domain().smallEnd(2);
2009 
2010  std::ofstream ofile;
2011  ofile.open ("MOST_xz_positions.txt");
2012 
2013  for (MFIter mfi(*fields[imf_cc], TileNoZ()); mfi.isValid(); ++mfi) {
2014  Box pbx = mfi.tilebox();
2015 
2016  if(pbx.smallEnd(2) != klo) { continue; }
2017 
2018  pbx.makeSlab(2,klo);
2019 
2020  int il = pbx.smallEnd(0); int iu = pbx.bigEnd(0);
2021 
2022  auto x_pos_arr = x_pos_mf->array(mfi);
2023  auto z_pos_arr = z_pos_mf->array(mfi);
2024 
2025  int k = 0;
2026  for (int i(il); i <= iu; ++i)
2027  ofile << x_pos_arr(i,j,k) << ' ' << z_pos_arr(i,j,k) << "\n";
2028  }
2029  ofile.close();
2030 }
Here is the call graph for this function:

Member Data Documentation

◆ include_subgrid_vel

bool MOSTAverage::include_subgrid_vel = false
protected

Referenced by make_MOSTAverage_at_level().

◆ m_averages

amrex::Vector<amrex::Vector<std::unique_ptr<amrex::MultiFab> > > MOSTAverage::m_averages
protected

Referenced by compute_eb_averages(), compute_plane_averages(), compute_region_averages(), get_average(), make_MOSTAverage_at_level(), set_plane_normalization(), and write_averages().

◆ m_eb_vec

amrex::Vector<const eb_*> MOSTAverage::m_eb_vec
protected

Referenced by compute_eb_averages(), and set_eb_normalization().

◆ m_fact_new

amrex::Real MOSTAverage::m_fact_new
protected

Referenced by compute_eb_averages(), compute_plane_averages(), compute_region_averages(), and make_MOSTAverage_at_level().

◆ m_fact_old

amrex::Real MOSTAverage::m_fact_old
protected

Referenced by compute_eb_averages(), compute_plane_averages(), compute_region_averages(), and make_MOSTAverage_at_level().

◆ m_fields

amrex::Vector<amrex::Vector<amrex::MultiFab*> > MOSTAverage::m_fields
protected

Referenced by compute_eb_averages(), compute_plane_averages(), compute_region_averages(), make_MOSTAverage_at_level(), set_k_indices_T(), set_norm_indices_T(), set_norm_positions_T(), set_plane_normalization(), set_rotated_fields(), set_z_positions_T(), update_field_ptrs(), write_averages(), write_k_indices(), write_norm_indices(), and write_xz_positions().

◆ m_geom

const amrex::Vector<amrex::Geometry> MOSTAverage::m_geom
protected

Referenced by compute_eb_averages(), compute_plane_averages(), compute_region_averages(), make_MOSTAverage_at_level(), set_eb_normalization(), set_k_indices_N(), set_k_indices_T(), set_norm_indices_T(), set_norm_positions_T(), set_plane_normalization(), set_rotated_fields(), set_z_positions_EB(), set_z_positions_T(), write_averages(), write_k_indices(), write_norm_indices(), and write_xz_positions().

◆ m_i_indx

amrex::Vector<std::unique_ptr<amrex::iMultiFab> > MOSTAverage::m_i_indx
protected

Referenced by compute_plane_averages(), compute_region_averages(), make_MOSTAverage_at_level(), set_norm_indices_T(), and write_norm_indices().

◆ m_interp

bool MOSTAverage::m_interp {false}
protected

Referenced by compute_plane_averages(), compute_region_averages(), and make_MOSTAverage_at_level().

◆ m_j_indx

amrex::Vector<std::unique_ptr<amrex::iMultiFab> > MOSTAverage::m_j_indx
protected

Referenced by compute_plane_averages(), compute_region_averages(), make_MOSTAverage_at_level(), set_norm_indices_T(), and write_norm_indices().

◆ m_k_in

amrex::Vector<int> MOSTAverage::m_k_in
protected

Referenced by set_k_indices_N(), and set_k_indices_T().

◆ m_k_indx

amrex::Vector<std::unique_ptr<amrex::iMultiFab> > MOSTAverage::m_k_indx
protected

Referenced by compute_plane_averages(), compute_region_averages(), get_k_indices(), make_MOSTAverage_at_level(), set_k_indices_N(), set_k_indices_T(), set_norm_indices_T(), write_k_indices(), and write_norm_indices().

◆ m_maxlev

int MOSTAverage::m_maxlev {0}
protected

Referenced by make_MOSTAverage_at_level(), set_eb_normalization(), and set_plane_normalization().

◆ m_mesh_type

MeshType MOSTAverage::m_mesh_type
protected

◆ m_navg

int MOSTAverage::m_navg {6}
protected

Referenced by compute_eb_averages(), compute_plane_averages(), compute_region_averages(), make_MOSTAverage_at_level(), set_eb_normalization(), set_plane_normalization(), and write_averages().

◆ m_ncell_plane

amrex::Vector<amrex::Vector<int> > MOSTAverage::m_ncell_plane
protected

Referenced by compute_plane_averages(), and set_plane_normalization().

◆ m_ncell_region

int MOSTAverage::m_ncell_region {1}
protected

Referenced by compute_region_averages(), and set_region_normalization().

◆ m_norm_vec

bool MOSTAverage::m_norm_vec {false}
protected

Referenced by make_MOSTAverage_at_level().

◆ m_nvar

int MOSTAverage::m_nvar {6}
protected

Referenced by make_MOSTAverage_at_level().

◆ m_plane_average

amrex::Vector<amrex::Vector<amrex::Real> > MOSTAverage::m_plane_average
protected

Referenced by compute_eb_averages(), compute_plane_averages(), set_eb_normalization(), and set_plane_normalization().

◆ m_policy

int MOSTAverage::m_policy {0}
protected

Referenced by compute_averages(), and make_MOSTAverage_at_level().

◆ m_pp_prefix

std::string MOSTAverage::m_pp_prefix
protected

Referenced by set_k_indices_N(), set_k_indices_T(), set_norm_indices_T(), set_norm_positions_T(), set_z_positions_EB(), and set_z_positions_T().

◆ m_radius

int MOSTAverage::m_radius {0}
protected

Referenced by compute_region_averages(), set_k_indices_N(), set_k_indices_T(), set_norm_indices_T(), and set_region_normalization().

◆ m_rot_fields

amrex::Vector<amrex::Vector<std::unique_ptr<amrex::MultiFab> > > MOSTAverage::m_rot_fields
protected

Referenced by compute_plane_averages(), compute_region_averages(), make_MOSTAverage_at_level(), and set_rotated_fields().

◆ m_rotate

bool MOSTAverage::m_rotate {false}
protected

Referenced by compute_averages(), compute_plane_averages(), compute_region_averages(), and make_MOSTAverage_at_level().

◆ m_t_avg

bool MOSTAverage::m_t_avg {false}
protected

Referenced by compute_averages(), compute_eb_averages(), compute_plane_averages(), compute_region_averages(), and make_MOSTAverage_at_level().

◆ m_t_init

amrex::Vector<int> MOSTAverage::m_t_init
protected

Referenced by compute_averages(), compute_eb_averages(), compute_plane_averages(), compute_region_averages(), and make_MOSTAverage_at_level().

◆ m_terrain_type

TerrainType MOSTAverage::m_terrain_type
protected

Referenced by make_MOSTAverage_at_level().

◆ m_time_window

amrex::Real MOSTAverage::m_time_window {amrex::Real(1.0e-16)}
protected

Referenced by make_MOSTAverage_at_level().

◆ m_total_bndry_area

amrex::Vector<amrex::Vector<amrex::Real> > MOSTAverage::m_total_bndry_area
protected

Referenced by compute_eb_averages(), and set_eb_normalization().

◆ m_Vsg

amrex::Vector<amrex::Real> MOSTAverage::m_Vsg
protected

Referenced by compute_eb_averages(), compute_plane_averages(), compute_region_averages(), and make_MOSTAverage_at_level().

◆ m_x_pos

amrex::Vector<std::unique_ptr<amrex::MultiFab> > MOSTAverage::m_x_pos
protected

Referenced by compute_plane_averages(), compute_region_averages(), make_MOSTAverage_at_level(), set_norm_positions_T(), set_z_positions_T(), and write_xz_positions().

◆ m_y_pos

amrex::Vector<std::unique_ptr<amrex::MultiFab> > MOSTAverage::m_y_pos
protected

Referenced by compute_plane_averages(), compute_region_averages(), make_MOSTAverage_at_level(), set_norm_positions_T(), and set_z_positions_T().

◆ m_z_phys_nd

amrex::Vector<amrex::MultiFab*> MOSTAverage::m_z_phys_nd
protected

Referenced by compute_plane_averages(), compute_region_averages(), make_MOSTAverage_at_level(), set_k_indices_T(), set_norm_indices_T(), set_norm_positions_T(), set_rotated_fields(), and set_z_positions_T().

◆ m_z_pos

amrex::Vector<std::unique_ptr<amrex::MultiFab> > MOSTAverage::m_z_pos
protected

Referenced by compute_plane_averages(), compute_region_averages(), make_MOSTAverage_at_level(), set_norm_positions_T(), set_z_positions_T(), and write_xz_positions().

◆ m_zref

amrex::Vector<std::unique_ptr<amrex::MultiFab> > MOSTAverage::m_zref
protected

Referenced by get_zref(), make_MOSTAverage_at_level(), set_k_indices_N(), set_k_indices_T(), set_norm_indices_T(), set_norm_positions_T(), set_z_positions_EB(), and set_z_positions_T().

◆ zref_default

const amrex::Real MOSTAverage::zref_default = amrex::Real(10.0)
protected

Referenced by make_MOSTAverage_at_level(), set_k_indices_N(), set_k_indices_T(), set_norm_indices_T(), set_norm_positions_T(), set_z_positions_EB(), and set_z_positions_T().

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