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::Orientation face, 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< amrex::Vector< amrex::Real >> &zlevels_stag, 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 validate_lateral_reference_stencil (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_average (const int &lev, const int &comp)
 
bool do_time_averaging () const
 
int get_navg () const
 
bool time_avg_is_initialized (const int &lev) const
 
void set_time_avg_initialized (const int &lev)
 
amrex::Vector< amrex::Real > get_plane_average (const int &lev) const
 
bool set_plane_average (const int &lev, const amrex::Vector< amrex::Real > &pavg)
 
amrex::MultiFab * get_zref (const int &lev) const
 
const amrex::iMultiFab * get_k_indices (const int &lev) const
 
amrex::IntVect get_ng_fill (const int &lev) const
 
void fill_planar_boundary (const int &lev, amrex::MultiFab &mf)
 
void extrap_ghost_cells (const int &lev, const int &iavg, const amrex::IntVect &ng_fill)
 

Static Public Member Functions

AMREX_GPU_HOST_DEVICE static AMREX_INLINE bool in_cell_z (const amrex::Real &z, const amrex::Real &z_lo, const amrex::Real &z_hi)
 
static int k_index_below (const amrex::Vector< amrex::Real > &zlevels, const amrex::Real &zref)
 
static amrex::Real cell_center_height (const amrex::Vector< amrex::Real > &zlevels, const int &k)
 
static amrex::Real column_height (const amrex::Vector< amrex::Real > &zlevels)
 
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

amrex::Orientation m_face
 
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
 
amrex::Vector< amrex::Vector< amrex::Real > > m_zlevels_stag
 
int m_nvar {6}
 
int m_navg {9}
 
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< PlanarBoundary > m_planar_bndry
 
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
 
double m_time_window {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)
 
const amrex::Real zref_sentinel = amrex::Real(-1.0)
 

Constructor & Destructor Documentation

◆ MOSTAverage() [1/3]

MOSTAverage::MOSTAverage ( amrex::Orientation  face,
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< amrex::Vector< amrex::Real >> &  zlevels_stag,
const amrex::Vector< const eb_ * > &  eb_vec = {} 
)
explicit

Construct the MOST averaging helper.

Parameters
[in]faceorientation of face for MOST averaging
[in]geomgeometry for all AMR levels
[in]has_zphyswhether nodal physical-height data are available
[in]a_pp_prefixParmParse prefix
[in]m_mesh_typemesh type
[in]m_terrain_typeterrain representation
[in]zlevels_stagnominal staggered z levels for all AMR levels
[in]eb_vecoptional embedded-boundary geometry data

◆ ~MOSTAverage()

MOSTAverage::~MOSTAverage ( )
inline

Destroy the MOST averaging helper.

45  {}

◆ MOSTAverage() [2/3]

MOSTAverage::MOSTAverage ( MOSTAverage &&  )
defaultnoexcept

Default move constructor.

◆ MOSTAverage() [3/3]

MOSTAverage::MOSTAverage ( const MOSTAverage &  other)
delete

Deleted copy constructor.

Parameters
[in]othersource object

Member Function Documentation

◆ cell_center_height()

static amrex::Real MOSTAverage::cell_center_height ( const amrex::Vector< amrex::Real > &  zlevels,
const int &  k 
)
inlinestatic

Height of a cell center above the bottom of the domain, on a column of staggered (face) heights.

Parameters
[in]zlevelsstaggered heights of the column (nz+1 faces)
[in]kcell index
386  {
387  return myhalf * (zlevels[k] + zlevels[k+1]) - zlevels[0];
388  }
myhalf
constexpr amrex::Real myhalf
Definition: ERF_NumericalConstants.H:41

Referenced by k_index_below(), and set_k_indices_N().

Here is the caller graph for this function:

◆ column_height()

static amrex::Real MOSTAverage::column_height ( const amrex::Vector< amrex::Real > &  zlevels)
inlinestatic

Height of the top of the column above its bottom, on a column of staggered (face) heights; the high z face's distances are measured from it.

Parameters
[in]zlevelsstaggered heights of the column (nz+1 faces)
397  {
398  return zlevels.back() - zlevels[0];
399  }

Referenced by set_k_indices_N().

Here is the caller graph for this function:

◆ compute_averages()

void MOSTAverage::compute_averages ( const int &  lev)

Driver for the selected average policy.

Parameters
[in]levlevel index

Function to call the type of average computation.

Parameters
[in]levCurrent level
1460 {
1461  if (m_rotate) set_rotated_fields(lev);
1462 
1463  switch(m_policy) {
1464  case 0: // Standard plane average
1465  compute_plane_averages(lev);
1466  break;
1467  case 1: // Local region/point
1468  compute_region_averages(lev);
1469  break;
1470  case 2: // EB Terrain average
1471  compute_eb_averages(lev);
1472  break;
1473  default:
1474  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(false, "Unknown policy for MOSTAverage!");
1475  }
1476 
1477  // We have initialized the averages
1478  if (m_t_avg) m_t_init[lev] = 1;
1479 }
AMREX_ALWAYS_ASSERT_WITH_MESSAGE
AMREX_ALWAYS_ASSERT_WITH_MESSAGE(axis >=0 &&axis< AMREX_SPACEDIM, "prob.axis must be 0, 1 or 2")
MOSTAverage::m_t_avg
bool m_t_avg
Definition: ERF_MOSTAverage.H:604
MOSTAverage::compute_plane_averages
void compute_plane_averages(const int &lev)
Definition: ERF_MOSTAverage.cpp:1488
MOSTAverage::m_policy
int m_policy
Definition: ERF_MOSTAverage.H:568
MOSTAverage::compute_region_averages
void compute_region_averages(const int &lev)
Definition: ERF_MOSTAverage.cpp:2087
MOSTAverage::m_t_init
amrex::Vector< int > m_t_init
Definition: ERF_MOSTAverage.H:605
MOSTAverage::set_rotated_fields
void set_rotated_fields(const int &lev)
Definition: ERF_MOSTAverage.cpp:513
MOSTAverage::m_rotate
bool m_rotate
Definition: ERF_MOSTAverage.H:569
MOSTAverage::compute_eb_averages
void compute_eb_averages(const int &lev)
Definition: ERF_MOSTAverage.cpp:2720
Here is the call graph for this function:

◆ compute_eb_averages()

void MOSTAverage::compute_eb_averages ( const int &  lev)

Fill averages for the embedded-boundary policy.

Parameters
[in]levlevel index

Function to compute average over an EB surface.

Parameters
[in]levCurrent level
2721 {
2722  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
2723  m_face.coordDir() == 2,
2724  "EB MOST averaging currently supports only z faces.");
2725  AMREX_ALWAYS_ASSERT(m_eb_vec[lev] != nullptr);
2726 
2727  // Peel back the level
2728  auto& fields = m_fields[lev];
2729  auto& averages = m_averages[lev];
2730  auto& plane_average = m_plane_average[lev];
2731 
2732  // Get EB data
2733  const auto& cc_flags = m_eb_vec[lev]->get_const_factory()->getMultiEBCellFlagFab();
2734  auto cc_afrac = m_eb_vec[lev]->get_const_factory()->getAreaFrac();
2735  const auto& cc_bnorm = m_eb_vec[lev]->get_const_factory()->getBndryNormal();
2736  const auto& u_vfrac = m_eb_vec[lev]->get_u_const_factory()->getVolFrac();
2737  const auto& v_vfrac = m_eb_vec[lev]->get_v_const_factory()->getVolFrac();
2738  const auto& w_vfrac = m_eb_vec[lev]->get_w_const_factory()->getVolFrac();
2739 
2740  // Get geometry for cell sizes
2741  auto const& dx_arr = m_geom[lev].CellSizeArray();
2742  Real dx = dx_arr[0];
2743  Real dy = dx_arr[1];
2744  Real dz = dx_arr[2];
2745 
2746  // Set factors for time averaging
2747  Real d_fact_new, d_fact_old;
2748  if (m_t_avg && m_t_init[lev]) {
2749  d_fact_new = m_fact_new;
2750  d_fact_old = m_fact_old;
2751  } else {
2752  d_fact_new = one;
2753  d_fact_old = zero;
2754  }
2755 
2756  // GPU array to accumulate averages into
2757  Gpu::DeviceVector<Real> pavg(plane_average.size(), zero);
2758  Real* plane_avg = pavg.data();
2759 
2760  // Vectors for normalization and buffer storage
2761  Vector<Real> denom(plane_average.size(),zero);
2762  Vector<Real> val_old(plane_average.size(),zero);
2763 
2764  //
2765  //----------------------------------------------------------
2766  // Averages for U, V, and tangential velocity (in local coordinate)
2767  //----------------------------------------------------------
2768  //
2769  {
2770  denom[0] = one / m_total_bndry_area[lev][0];
2771  val_old[0] = plane_average[0]*d_fact_old;
2772  denom[1] = one / m_total_bndry_area[lev][1];
2773  val_old[1] = plane_average[1]*d_fact_old;
2774 
2775  int iavg = m_navg - 1; // Tangential velocity magnitude
2776  denom[iavg] = one / m_total_bndry_area[lev][iavg];
2777  val_old[iavg] = plane_average[iavg]*d_fact_old;
2778 
2779  const Real Vsg = m_Vsg[lev]; // Subgrid scale velocity
2780 
2781  //
2782  // NOTE: deliberately not threaded. The body is a pure reduction into the shared
2783  // plane_avg entries, and on the host Gpu::deviceReduceSum is an unordered
2784  // "#pragma omp atomic" add. Threading it would make the surface-layer averages
2785  // -- and through them the surface fluxes -- depend on thread scheduling in the
2786  // last bits. See the same note in PlaneAverage::compute_averages.
2787  //
2788  for (MFIter mfi(*fields[3], TileNoZ()); mfi.isValid(); ++mfi) {
2789  const auto& flag = cc_flags[mfi];
2790 
2791  // Skip boxes that are not singlevalued (MultiCutFab only has data for singlevalued boxes)
2792  if (flag.getType() != FabType::singlevalued) continue;
2793 
2794  Box bx = mfi.tilebox(); // Full 3D box
2795 
2796  // Get EB arrays
2797  auto const flag_arr = flag.const_array();
2798  auto const afrac_x = cc_afrac[0]->const_array(mfi);
2799  auto const afrac_y = cc_afrac[1]->const_array(mfi);
2800  auto const afrac_z = cc_afrac[2]->const_array(mfi);
2801  auto const bnorm_arr = cc_bnorm.const_array(mfi);
2802  auto const u_vf_arr = u_vfrac.const_array(mfi);
2803  auto const v_vf_arr = v_vfrac.const_array(mfi);
2804  auto const w_vf_arr = w_vfrac.const_array(mfi);
2805 
2806  // Get velocity arrays
2807  auto const u_arr = fields[0]->const_array(mfi);
2808  auto const v_arr = fields[1]->const_array(mfi);
2809  auto const w_arr = fields[2]->const_array(mfi);
2810 
2811  ParallelFor(Gpu::KernelInfo().setReduction(true), bx, [=]
2812  AMREX_GPU_DEVICE(int i, int j, int k, Gpu::Handler const& handler) noexcept
2813  {
2814  // Area-weighted averaging over cut cells at any k
2815  if (flag_arr(i,j,k).isSingleValued()) {
2816  // Compute area from face-centered area fractions
2817  Real axm = afrac_x(i ,j ,k );
2818  Real axp = afrac_x(i+1,j ,k );
2819  Real aym = afrac_y(i ,j ,k );
2820  Real ayp = afrac_y(i ,j+1,k );
2821  Real azm = afrac_z(i ,j ,k );
2822  Real azp = afrac_z(i ,j ,k+1);
2823 
2824  Real adx = (axm - axp) * dy * dz;
2825  Real ady = (aym - ayp) * dx * dz;
2826  Real adz = (azm - azp) * dx * dy;
2827 
2828  Real area = std::sqrt(adx*adx + ady*ady + adz*adz);
2829 
2830  // Volume-weighted interpolation of velocities to cell center
2831  Real vf_u_lo = u_vf_arr(i,j,k);
2832  Real vf_u_hi = u_vf_arr(i+1,j,k);
2833  Real sum_vf_u = vf_u_lo + vf_u_hi;
2834  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;
2835 
2836  Real vf_v_lo = v_vf_arr(i,j,k);
2837  Real vf_v_hi = v_vf_arr(i,j+1,k);
2838  Real sum_vf_v = vf_v_lo + vf_v_hi;
2839  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;
2840 
2841  Real vf_w_lo = w_vf_arr(i,j,k);
2842  Real vf_w_hi = w_vf_arr(i,j,k+1);
2843  Real sum_vf_w = vf_w_lo + vf_w_hi;
2844  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;
2845 
2846  // Get normal vector components
2847  Real nx = bnorm_arr(i,j,k,0);
2848  Real ny = bnorm_arr(i,j,k,1);
2849  Real nz = bnorm_arr(i,j,k,2);
2850 
2851  // Compute tangential velocity components
2852  Real v_dot_n = u_cc*nx + v_cc*ny + w_cc*nz;
2853  Real u_tangent = u_cc - v_dot_n * nx;
2854  Real v_tangent = v_cc - v_dot_n * ny;
2855  Real mag = std::sqrt(u_tangent*u_tangent + v_tangent*v_tangent + Vsg*Vsg);
2856 
2857  // Area-weighted sum
2858  Real val_u = u_tangent * area;
2859  Real val_v = v_tangent * area;
2860  Real val_mag = mag * area;
2861 
2862  Gpu::deviceReduceSum(&plane_avg[0], val_u, handler);
2863  Gpu::deviceReduceSum(&plane_avg[1], val_v, handler);
2864  Gpu::deviceReduceSum(&plane_avg[iavg], val_mag, handler);
2865  }
2866  });
2867  }
2868  }
2869 
2870  //
2871  //----------------------------------------------------------
2872  // Averages for T,Qv (cell-centered scalars)
2873  //----------------------------------------------------------
2874  //
2875  for (int imf(3); imf < 5; ++imf) {
2876 
2877  // Continue if no valid Qv pointer
2878  if (!fields[imf]) continue;
2879 
2880  denom[imf] = one / m_total_bndry_area[lev][imf];
2881  val_old[imf] = plane_average[imf]*d_fact_old;
2882 
2883  //
2884  // NOTE: deliberately not threaded. The body is a pure reduction into the shared
2885  // plane_avg entries, and on the host Gpu::deviceReduceSum is an unordered
2886  // "#pragma omp atomic" add. Threading it would make the surface-layer averages
2887  // -- and through them the surface fluxes -- depend on thread scheduling in the
2888  // last bits. See the same note in PlaneAverage::compute_averages.
2889  //
2890  for (MFIter mfi(*fields[imf], TileNoZ()); mfi.isValid(); ++mfi) {
2891  const auto& flag = cc_flags[mfi];
2892 
2893  // Skip boxes that are not singlevalued (MultiCutFab only has data for singlevalued boxes)
2894  if (flag.getType() != FabType::singlevalued) continue;
2895 
2896  Box bx = mfi.tilebox(); // Full 3D box
2897  auto const flag_arr = flag.const_array();
2898  auto const afrac_x = cc_afrac[0]->const_array(mfi);
2899  auto const afrac_y = cc_afrac[1]->const_array(mfi);
2900  auto const afrac_z = cc_afrac[2]->const_array(mfi);
2901  auto const mf_arr = fields[imf]->const_array(mfi);
2902 
2903  ParallelFor(Gpu::KernelInfo().setReduction(true), bx, [=]
2904  AMREX_GPU_DEVICE(int i, int j, int k, Gpu::Handler const& handler) noexcept
2905  {
2906  // Area-weighted averaging over cut cells at any k
2907  if (flag_arr(i,j,k).isSingleValued()) {
2908  // Compute area from face-centered area fractions
2909  Real axm = afrac_x(i ,j ,k );
2910  Real axp = afrac_x(i+1,j ,k );
2911  Real aym = afrac_y(i ,j ,k );
2912  Real ayp = afrac_y(i ,j+1,k );
2913  Real azm = afrac_z(i ,j ,k );
2914  Real azp = afrac_z(i ,j ,k+1);
2915 
2916  Real adx = (axm - axp) * dy * dz;
2917  Real ady = (aym - ayp) * dx * dz;
2918  Real adz = (azm - azp) * dx * dy;
2919 
2920  Real area = std::sqrt(adx*adx + ady*ady + adz*adz);
2921 
2922  Real val = mf_arr(i,j,k) * area;
2923  Gpu::deviceReduceSum(&plane_avg[imf], val, handler);
2924  }
2925  });
2926  }
2927  }
2928 
2929  //
2930  //------------------------------------------------------------------------
2931  // Averages for virtual potential temperature
2932  //------------------------------------------------------------------------
2933  //
2934  if (fields[4]) // We have water vapor
2935  {
2936  int iavg = 5;
2937  denom[iavg] = one / m_total_bndry_area[lev][iavg];
2938  val_old[iavg] = plane_average[iavg]*d_fact_old;
2939 
2940  //
2941  // NOTE: deliberately not threaded. The body is a pure reduction into the shared
2942  // plane_avg entries, and on the host Gpu::deviceReduceSum is an unordered
2943  // "#pragma omp atomic" add. Threading it would make the surface-layer averages
2944  // -- and through them the surface fluxes -- depend on thread scheduling in the
2945  // last bits. See the same note in PlaneAverage::compute_averages.
2946  //
2947  for (MFIter mfi(*fields[4], TileNoZ()); mfi.isValid(); ++mfi)
2948  {
2949  const auto& flag = cc_flags[mfi];
2950 
2951  // Skip boxes that are not singlevalued (MultiCutFab only has data for singlevalued boxes)
2952  if (flag.getType() != FabType::singlevalued) continue;
2953 
2954  Box bx = mfi.tilebox(); // Full 3D box
2955  auto const flag_arr = flag.const_array();
2956  auto const afrac_x = cc_afrac[0]->const_array(mfi);
2957  auto const afrac_y = cc_afrac[1]->const_array(mfi);
2958  auto const afrac_z = cc_afrac[2]->const_array(mfi);
2959 
2960  const Array4<Real const> T_mf_arr = fields[3]->const_array(mfi);
2961  const Array4<Real const> qv_mf_arr = fields[4]->const_array(mfi);
2962  const Array4<Real const> qr_mf_arr = (fields[5]) ? fields[5]->const_array(mfi) :
2963  Array4<const Real> {};
2964 
2965  ParallelFor(Gpu::KernelInfo().setReduction(true), bx, [=]
2966  AMREX_GPU_DEVICE(int i, int j, int k, Gpu::Handler const& handler) noexcept
2967  {
2968  if (flag_arr(i,j,k).isSingleValued()) {
2969  // Compute area from face-centered area fractions
2970  Real axm = afrac_x(i ,j ,k );
2971  Real axp = afrac_x(i+1,j ,k );
2972  Real aym = afrac_y(i ,j ,k );
2973  Real ayp = afrac_y(i ,j+1,k );
2974  Real azm = afrac_z(i ,j ,k );
2975  Real azp = afrac_z(i ,j ,k+1);
2976 
2977  Real adx = (axm - axp) * dy * dz;
2978  Real ady = (aym - ayp) * dx * dz;
2979  Real adz = (azm - azp) * dx * dy;
2980 
2981  Real area = std::sqrt(adx*adx + ady*ady + adz*adz);
2982 
2983  Real vfac;
2984  if (qr_mf_arr) {
2985  // We also have liquid water
2986  vfac = one + epsv*qv_mf_arr(i,j,k) - qr_mf_arr(i,j,k);
2987  } else {
2988  vfac = one + epsv*qv_mf_arr(i,j,k);
2989  }
2990  const Real val = T_mf_arr(i,j,k) * vfac * area;
2991  Gpu::deviceReduceSum(&plane_avg[iavg], val, handler);
2992  }
2993  });
2994  }
2995  }
2996  else // copy temperature
2997  {
2998  int iavg = m_navg - 4;
2999  denom[iavg] = one / m_total_bndry_area[lev][iavg];
3000  // plane_avg[iavg] = plane_avg[2]
3001  Gpu::copy(Gpu::deviceToDevice, pavg.begin() + 3, pavg.begin() + 4,
3002  pavg.begin() + iavg);
3003  }
3004 
3005  // Copy to host and sum across procs
3006  Gpu::copy(Gpu::deviceToHost, pavg.begin(), pavg.end(), plane_average.begin());
3007  ParallelDescriptor::ReduceRealSum(plane_average.data(), plane_average.size());
3008 
3009  // Normalize by total area and apply time averaging
3010  for (int iavg(0); iavg < m_navg; ++iavg){
3011  plane_average[iavg] *= denom[iavg]*d_fact_new;
3012  plane_average[iavg] += val_old[iavg];
3013  averages[iavg]->setVal(plane_average[iavg]);
3014  }
3015 }
epsv
constexpr amrex::Real epsv
Definition: ERF_Constants.H:56
vfac
Real vfac
Definition: ERF_InitCustomPertVels_ABL.H:42
nx
const int nx
Definition: ERF_InitCustomPertVels_CloudChamber.H:14
ny
const int ny
Definition: ERF_InitCustomPertVels_CloudChamber.H:15
dy
const Real dy
Definition: ERF_InitCustomPert_ABL.H:45
dx
const Real dx
Definition: ERF_InitCustomPert_ABL.H:44
AMREX_ALWAYS_ASSERT
AMREX_ALWAYS_ASSERT(bx.length()[2]==khi+1)
ParallelFor
ParallelFor(fab_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);})
one
constexpr amrex::Real one
Definition: ERF_NumericalConstants.H:37
zero
constexpr amrex::Real zero
Definition: ERF_NumericalConstants.H:36
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:566
MOSTAverage::m_averages
amrex::Vector< amrex::Vector< std::unique_ptr< amrex::MultiFab > > > m_averages
Definition: ERF_MOSTAverage.H:577
MOSTAverage::m_total_bndry_area
amrex::Vector< amrex::Vector< amrex::Real > > m_total_bndry_area
Definition: ERF_MOSTAverage.H:600
MOSTAverage::m_face
amrex::Orientation m_face
Definition: ERF_MOSTAverage.H:554
MOSTAverage::m_Vsg
amrex::Vector< amrex::Real > m_Vsg
Definition: ERF_MOSTAverage.H:612
MOSTAverage::m_plane_average
amrex::Vector< amrex::Vector< amrex::Real > > m_plane_average
Definition: ERF_MOSTAverage.H:584
MOSTAverage::m_fact_new
amrex::Real m_fact_new
Definition: ERF_MOSTAverage.H:607
MOSTAverage::m_fields
amrex::Vector< amrex::Vector< amrex::MultiFab * > > m_fields
Definition: ERF_MOSTAverage.H:556
MOSTAverage::m_eb_vec
amrex::Vector< const eb_ * > m_eb_vec
Definition: ERF_MOSTAverage.H:599
MOSTAverage::m_fact_old
amrex::Real m_fact_old
Definition: ERF_MOSTAverage.H:607
MOSTAverage::m_geom
const amrex::Vector< amrex::Geometry > m_geom
Definition: ERF_MOSTAverage.H:555
WDM6SedCellScratch::dz
@ dz
Definition: ERF_AdvanceWDM6.cpp:273

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)

Fill averages for the plane policy.

Parameters
[in]levlevel index

Function to compute average over a plane.

Parameters
[in]levCurrent level
1489 {
1490  // Peel back the level
1491  auto& fields = m_fields[lev];
1492  auto& rot_fields = m_rot_fields[lev];
1493  auto& averages = m_averages[lev];
1494  const auto & geom = m_geom[lev];
1495 
1496  auto& z_phys = m_z_phys_nd[lev];
1497  auto& x_pos = m_x_pos[lev];
1498  auto& y_pos = m_y_pos[lev];
1499  auto& z_pos = m_z_pos[lev];
1500 
1501  auto& i_indx = m_i_indx[lev];
1502  auto& j_indx = m_j_indx[lev];
1503  auto& k_indx = m_k_indx[lev];
1504 
1505  auto& ncell_plane = m_ncell_plane[lev];
1506  auto& plane_average = m_plane_average[lev];
1507 
1508  // Set factors for time averaging
1509  Real d_fact_new, d_fact_old;
1510  if (m_t_avg && m_t_init[lev]) {
1511  d_fact_new = m_fact_new;
1512  d_fact_old = m_fact_old;
1513  } else {
1514  d_fact_new = one;
1515  d_fact_old = zero;
1516  }
1517 
1518 
1519  // GPU array to accumulate averages into
1520  Gpu::DeviceVector<Real> pavg(plane_average.size(), zero);
1521  Real* plane_avg = pavg.data();
1522 
1523  // Vectors for normalization and buffer storage
1524  Vector<Real> denom(plane_average.size(),zero);
1525  Vector<Real> val_old(plane_average.size(),zero);
1526 
1527  //
1528  //----------------------------------------------------------
1529  // Averages over all the fields
1530  //----------------------------------------------------------
1531  //
1532  Box domain = geom.Domain();
1533 
1534  Array<int,AMREX_SPACEDIM> is_per = {0,0,0};
1535  for (int idim(0); idim < AMREX_SPACEDIM; ++idim) {
1536  if (geom.isPeriodic(idim)) is_per[idim] = 1;
1537  }
1538 
1539  const int dir = m_face.coordDir();
1540  // The non-terrain reference is constant on a level. Do not index the
1541  // cell-centered k-index FAB with an MFIter from a staggered velocity FAB;
1542  // distributed BoxArrays have the same box ordering but different valid
1543  // extents. Spatial index arrays are retained for terrain-normal cases.
1544  const bool fitted_terrain =
1545  (m_terrain_type == TerrainType::StaticFittedMesh) ||
1546  (m_terrain_type == TerrainType::MovingFittedMesh);
1547  const bool use_spatial_indices =
1548  m_norm_vec || (dir == 2 && fitted_terrain);
1549  if (!use_spatial_indices) {
1550  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
1551  m_k_indx[lev] != nullptr,
1552  "Plane averaging requires a reference-index field.");
1553  }
1554  const int wall_normal_ref = use_spatial_indices ? 0 : m_k_indx[lev]->min(0);
1555  // Averages for U,V,W,T,Qv (not Qc)
1556  for (int imf(0); imf < 5; ++imf) {
1557 
1558  int sm_index = 0;
1559  if (m_face.isLow()) {
1560  sm_index = m_geom[lev].Domain().smallEnd(dir);
1561  } else {
1562  sm_index = m_geom[lev].Domain().bigEnd(dir);
1563  if (imf < 3 && imf == dir) {
1564  sm_index += 1;
1565  }
1566  }
1567 
1568  const int normal_face_offset =
1569  (!m_face.isLow() && imf < 3 && imf == dir) ? 1 : 0;
1570 
1571  // Continue if no valid Qv pointer
1572  if (!fields[imf]) continue;
1573 
1574  denom[imf] = one / (Real)ncell_plane[imf];
1575  val_old[imf] = plane_average[imf]*d_fact_old;
1576 
1577  //
1578  // NOTE: deliberately not threaded. The body is a pure reduction into the shared
1579  // plane_avg entries, and on the host Gpu::deviceReduceSum is an unordered
1580  // "#pragma omp atomic" add. Threading it would make the surface-layer averages
1581  // -- and through them the surface fluxes -- depend on thread scheduling in the
1582  // last bits. See the same note in PlaneAverage::compute_averages.
1583  //
1584  for (MFIter mfi(*fields[imf], TileNoZ()); mfi.isValid(); ++mfi) {
1585  Box vbx = mfi.validbox(); // This is the grid (not tile)
1586  Box pbx = mfi.tilebox(); // This is the tile (not grid)
1587 
1588  if (use_spatial_indices) {
1589  if (m_face.isLow()) {
1590  if (vbx.smallEnd(dir) != sm_index ||
1591  pbx.smallEnd(dir) != sm_index) {
1592  continue;
1593  }
1594  } else {
1595  if (vbx.bigEnd(dir) != sm_index ||
1596  pbx.bigEnd(dir) != sm_index) {
1597  continue;
1598  }
1599  }
1600  } else {
1601  // For Cartesian averages, the source FAB must own the
1602  // requested reference plane, not necessarily the boundary
1603  // face. A distributed high-side FAB can own the face while
1604  // the reference plane is in another FAB.
1605  Box reference_box = pbx;
1606  reference_box.setRange(dir, wall_normal_ref);
1607  if (!pbx.contains(reference_box)) {
1608  continue;
1609  }
1610  }
1611 
1612  // Make planar since mfiter is over fields
1613  pbx.setRange(dir, use_spatial_indices ? sm_index : wall_normal_ref);
1614 
1615  // Avoid double counting nodal data by changing the high end when we are
1616  // at the high side of the grid (not just of the tile)
1617  IndexType ixt = averages[imf]->boxArray().ixType();
1618  for (int idim(0); idim < AMREX_SPACEDIM; ++idim) {
1619  if (idim != dir && ixt.nodeCentered(idim) &&
1620  pbx.bigEnd(idim) == vbx.bigEnd(idim)) {
1621  int dom_hi = domain.bigEnd(idim) + 1;
1622  if (pbx.bigEnd(idim) < dom_hi || is_per[idim]) {
1623  pbx.growHi(idim,-1);
1624  }
1625  }
1626  }
1627 
1628  auto mf_arr = (m_rotate && imf != 2) ? rot_fields[imf]->const_array(mfi) :
1629  fields[imf]->const_array(mfi);
1630 
1631  if (m_interp) {
1632  const auto plo = geom.ProbLoArray();
1633  const auto dxInv = geom.InvCellSizeArray();
1634  const auto z_phys_arr = z_phys->const_array(mfi);
1635  auto x_pos_arr = x_pos->array(mfi);
1636  auto y_pos_arr = y_pos->array(mfi);
1637  auto z_pos_arr = z_pos->array(mfi);
1638  ParallelFor(Gpu::KernelInfo().setReduction(true), pbx, [=]
1639  AMREX_GPU_DEVICE(int i, int j, int k, Gpu::Handler const& handler) noexcept
1640  {
1641  Real interp{0};
1642  trilinear_interp_T(x_pos_arr(i,j,k), y_pos_arr(i,j,k), z_pos_arr(i,j,k),
1643  &interp, mf_arr, z_phys_arr, plo, dxInv, 1);
1644  Real val = interp;
1645  Gpu::deviceReduceSum(&plane_avg[imf], val, handler);
1646  });
1647  } else {
1648  auto k_arr = use_spatial_indices
1649  ? k_indx->const_array(mfi) : Array4<const int>{};
1650  auto j_arr = j_indx ? j_indx->const_array(mfi) : Array4<const int> {};
1651  auto i_arr = i_indx ? i_indx->const_array(mfi) : Array4<const int> {};
1652  // Spatial index fields are populated in their ghost cells
1653  // for local averaging. Clamp to the FAB extent, not just
1654  // the valid box, so those terrain-dependent ghost indices
1655  // remain available.
1656  const Box k_box = use_spatial_indices
1657  ? k_indx->fabbox(mfi.index()) : Box{};
1658  ParallelFor(Gpu::KernelInfo().setReduction(true), pbx, [=]
1659  AMREX_GPU_DEVICE(int i, int j, int k, Gpu::Handler const& handler) noexcept
1660  {
1661  const int ki = use_spatial_indices
1662  ? max(k_box.smallEnd(0), min(k_box.bigEnd(0), i)) : i;
1663  const int kj = use_spatial_indices
1664  ? max(k_box.smallEnd(1), min(k_box.bigEnd(1), j)) : j;
1665  const int kk = use_spatial_indices
1666  ? max(k_box.smallEnd(2), min(k_box.bigEnd(2), k)) : k;
1667  const int ref = (use_spatial_indices
1668  ? k_arr(ki,kj,kk) : wall_normal_ref) + normal_face_offset;
1669  int mi = i_arr ? i_arr(ki,kj,k) : i;
1670  int mj = j_arr ? j_arr(ki,kj,k) : j;
1671  int mk = k;
1672  if (dir == 0) {
1673  mi = ref;
1674  } else if (dir == 1) {
1675  mj = ref;
1676  } else {
1677  mk = ref;
1678  }
1679  Real val = mf_arr(mi,mj,mk);
1680  Gpu::deviceReduceSum(&plane_avg[imf], val, handler);
1681  });
1682  }
1683  }
1684  }
1685 
1686  //
1687  //------------------------------------------------------------------------
1688  // Averages for virtual potential temperature
1689  // (This is cell-centered so we don't need to worry about double-counting)
1690  //------------------------------------------------------------------------
1691  //
1692  if (fields[4]) // We have water vapor
1693  {
1694  int sm_index = 0;
1695  if (m_face.isLow()) {
1696  sm_index = m_geom[lev].Domain().smallEnd(dir);
1697  } else {
1698  sm_index = m_geom[lev].Domain().bigEnd(dir);
1699  }
1700 
1701  int iavg = 5;
1702  denom[iavg] = one / (Real)ncell_plane[iavg];
1703  val_old[iavg] = plane_average[iavg]*d_fact_old;
1704 
1705  //
1706  // NOTE: deliberately not threaded. The body is a pure reduction into the shared
1707  // plane_avg entries, and on the host Gpu::deviceReduceSum is an unordered
1708  // "#pragma omp atomic" add. Threading it would make the surface-layer averages
1709  // -- and through them the surface fluxes -- depend on thread scheduling in the
1710  // last bits. See the same note in PlaneAverage::compute_averages.
1711  //
1712  for (MFIter mfi(*fields[4], TileNoZ()); mfi.isValid(); ++mfi)
1713  {
1714  Box pbx = mfi.tilebox();
1715  const Box vbx = mfi.validbox();
1716 
1717  if (use_spatial_indices) {
1718  if (m_face.isLow()) {
1719  if (vbx.smallEnd(dir) != sm_index ||
1720  pbx.smallEnd(dir) != sm_index) {
1721  continue;
1722  }
1723  } else {
1724  if (vbx.bigEnd(dir) != sm_index ||
1725  pbx.bigEnd(dir) != sm_index) {
1726  continue;
1727  }
1728  }
1729  } else {
1730  Box reference_box = pbx;
1731  reference_box.setRange(dir, wall_normal_ref);
1732  if (!pbx.contains(reference_box)) {
1733  continue;
1734  }
1735  }
1736 
1737  pbx.setRange(dir, use_spatial_indices ? sm_index : wall_normal_ref);
1738 
1739  const Array4<Real const>& T_mf_arr = fields[3]->const_array(mfi);
1740  const Array4<Real const>& qv_mf_arr = (fields[4])? fields[4]->const_array(mfi) : Array4<const Real>{};
1741  const Array4<Real const>& qr_mf_arr = (fields[5])? fields[5]->const_array(mfi) : Array4<const Real>{};
1742 
1743  if (m_interp) {
1744  const auto plo = m_geom[lev].ProbLoArray();
1745  const auto dxInv = m_geom[lev].InvCellSizeArray();
1746  const auto z_phys_arr = z_phys->const_array(mfi);
1747  auto x_pos_arr = x_pos->array(mfi);
1748  auto y_pos_arr = y_pos->array(mfi);
1749  auto z_pos_arr = z_pos->array(mfi);
1750  ParallelFor(Gpu::KernelInfo().setReduction(true), pbx, [=]
1751  AMREX_GPU_DEVICE(int i, int j, int k, Gpu::Handler const& handler) noexcept
1752  {
1753  Real T_interp{0};
1754  Real qv_interp{0};
1755  trilinear_interp_T(x_pos_arr(i,j,k), y_pos_arr(i,j,k), z_pos_arr(i,j,k),
1756  &T_interp, T_mf_arr, z_phys_arr, plo, dxInv, 1);
1757  trilinear_interp_T(x_pos_arr(i,j,k), y_pos_arr(i,j,k), z_pos_arr(i,j,k),
1758  &qv_interp, qv_mf_arr, z_phys_arr, plo, dxInv, 1);
1759  Real vfac;
1760  if (qr_mf_arr) {
1761  // We also have liquid water
1762  Real qr_interp{0};
1763  trilinear_interp_T(x_pos_arr(i,j,k), y_pos_arr(i,j,k), z_pos_arr(i,j,k),
1764  &qr_interp, qr_mf_arr, z_phys_arr, plo, dxInv, 1);
1765  vfac = one + epsv*qv_interp - qr_interp;
1766  } else {
1767  vfac = one + epsv*qv_interp;
1768  }
1769  const Real val = T_interp * vfac;
1770  Gpu::deviceReduceSum(&plane_avg[iavg], val, handler);
1771  });
1772  } else {
1773  auto k_arr = use_spatial_indices
1774  ? k_indx->const_array(mfi) : Array4<const int>{};
1775  auto j_arr = j_indx ? j_indx->const_array(mfi) : Array4<const int> {};
1776  auto i_arr = i_indx ? i_indx->const_array(mfi) : Array4<const int> {};
1777  const Box k_box = use_spatial_indices
1778  ? k_indx->fabbox(mfi.index()) : Box{};
1779  ParallelFor(Gpu::KernelInfo().setReduction(true), pbx, [=]
1780  AMREX_GPU_DEVICE(int i, int j, int k, Gpu::Handler const& handler) noexcept
1781  {
1782  const int ki = use_spatial_indices
1783  ? max(k_box.smallEnd(0), min(k_box.bigEnd(0), i)) : i;
1784  const int kj = use_spatial_indices
1785  ? max(k_box.smallEnd(1), min(k_box.bigEnd(1), j)) : j;
1786  const int kk = use_spatial_indices
1787  ? max(k_box.smallEnd(2), min(k_box.bigEnd(2), k)) : k;
1788  const int ref = use_spatial_indices
1789  ? k_arr(ki,kj,kk) : wall_normal_ref;
1790  int mi = i_arr ? i_arr(ki,kj,k) : i;
1791  int mj = j_arr ? j_arr(ki,kj,k) : j;
1792  int mk = k;
1793  if (dir == 0) {
1794  mi = ref;
1795  } else if (dir == 1) {
1796  mj = ref;
1797  } else {
1798  mk = ref;
1799  }
1800  Real vfac;
1801  if (qr_mf_arr) {
1802  // We also have liquid water
1803  vfac = one + epsv*qv_mf_arr(mi,mj,mk) - qr_mf_arr(mi,mj,mk);
1804  } else {
1805  vfac = one + epsv*qv_mf_arr(mi,mj,mk);
1806  }
1807  const Real val = T_mf_arr(mi,mj,mk) * vfac;
1808  Gpu::deviceReduceSum(&plane_avg[iavg], val, handler);
1809  });
1810  }
1811  }
1812  }
1813  else // copy temperature
1814  {
1815  int iavg = m_navg - 4;
1816  denom[iavg] = one / (Real)ncell_plane[iavg];
1817  // plane_avg[iavg] = plane_avg[2]
1818  Gpu::copy(Gpu::deviceToDevice, pavg.begin() + 3, pavg.begin() + 4,
1819  pavg.begin() + iavg);
1820  }
1821 
1822  //
1823  //------------------------------------------------------------------------
1824  // Averages for the tangential velocity magnitude
1825  // (This is cell-centered so we don't need to worry about double-counting)
1826  //------------------------------------------------------------------------
1827  //
1828  {
1829  int sm_index = 0;
1830  if (m_face.isLow()) {
1831  sm_index = m_geom[lev].Domain().smallEnd(dir);
1832  } else {
1833  sm_index = m_geom[lev].Domain().bigEnd(dir);
1834  }
1835 
1836  const int imf_cc = 3;
1837  const int imf = 0;
1838  const int iavg = m_navg - 3;
1839  const int iavg_xz = m_navg - 2;
1840  const int iavg_yz = m_navg - 1;
1841  denom[iavg] = one / (Real)ncell_plane[iavg];
1842  denom[iavg_xz] = one / (Real)ncell_plane[iavg_xz];
1843  denom[iavg_yz] = one / (Real)ncell_plane[iavg_yz];
1844  val_old[iavg] = plane_average[iavg]*d_fact_old;
1845  val_old[iavg_xz] = plane_average[iavg_xz]*d_fact_old;
1846  val_old[iavg_yz] = plane_average[iavg_yz]*d_fact_old;
1847 
1848  const Real Vsg = m_Vsg[lev];
1849 
1850  //
1851  // NOTE: deliberately not threaded. The body is a pure reduction into the shared
1852  // plane_avg entries, and on the host Gpu::deviceReduceSum is an unordered
1853  // "#pragma omp atomic" add. Threading it would make the surface-layer averages
1854  // -- and through them the surface fluxes -- depend on thread scheduling in the
1855  // last bits. See the same note in PlaneAverage::compute_averages.
1856  //
1857  for (MFIter mfi(*fields[imf_cc], TileNoZ()); mfi.isValid(); ++mfi)
1858  {
1859  Box pbx = mfi.tilebox();
1860  const Box vbx = mfi.validbox();
1861 
1862  if (use_spatial_indices) {
1863  if (m_face.isLow()) {
1864  if (vbx.smallEnd(dir) != sm_index ||
1865  pbx.smallEnd(dir) != sm_index) {
1866  continue;
1867  }
1868  } else {
1869  if (vbx.bigEnd(dir) != sm_index ||
1870  pbx.bigEnd(dir) != sm_index) {
1871  continue;
1872  }
1873  }
1874  } else {
1875  Box reference_box = pbx;
1876  reference_box.setRange(dir, wall_normal_ref);
1877  if (!pbx.contains(reference_box)) {
1878  continue;
1879  }
1880  }
1881 
1882  pbx.setRange(dir, use_spatial_indices ? sm_index : wall_normal_ref);
1883 
1884  // Last element is Umag and always cell centered
1885  auto u_mf_arr = (m_rotate) ? rot_fields[imf ]->const_array(mfi) :
1886  fields[imf ]->const_array(mfi);
1887  auto v_mf_arr = (m_rotate) ? rot_fields[imf+1]->const_array(mfi) :
1888  fields[imf+1]->const_array(mfi);
1889  auto w_mf_arr = fields[imf+2]->const_array(mfi);
1890 
1891  if (m_interp) {
1892  // TODO: trilinear interp needs to be fixed for X and Y faces
1893  const auto plo = m_geom[lev].ProbLoArray();
1894  const auto dxInv = m_geom[lev].InvCellSizeArray();
1895  const auto z_phys_arr = z_phys->const_array(mfi);
1896  auto x_pos_arr = x_pos->array(mfi);
1897  auto y_pos_arr = y_pos->array(mfi);
1898  auto z_pos_arr = z_pos->array(mfi);
1899  ParallelFor(Gpu::KernelInfo().setReduction(true), pbx, [=]
1900  AMREX_GPU_DEVICE(int i, int j, int k, Gpu::Handler const& handler) noexcept
1901  {
1902  Real u_interp{0};
1903  Real v_interp{0};
1904  trilinear_interp_T(x_pos_arr(i,j,k), y_pos_arr(i,j,k), z_pos_arr(i,j,k),
1905  &u_interp, u_mf_arr, z_phys_arr, plo, dxInv, 1);
1906  trilinear_interp_T(x_pos_arr(i,j,k), y_pos_arr(i,j,k), z_pos_arr(i,j,k),
1907  &v_interp, v_mf_arr, z_phys_arr, plo, dxInv, 1);
1908  const Real val = std::sqrt(u_interp*u_interp + v_interp*v_interp + Vsg*Vsg);
1909  Gpu::deviceReduceSum(&plane_avg[iavg], val, handler);
1910 
1911  if (dir < 2) {
1912  // averages for mean velocity on XZ and YZ planes
1913  Real w_interp{0};
1914  trilinear_interp_T(x_pos_arr(i,j,k), y_pos_arr(i,j,k), z_pos_arr(i,j,k),
1915  &w_interp, w_mf_arr, z_phys_arr, plo, dxInv, 1);
1916  const Real val_xz = std::sqrt(u_interp*u_interp + w_interp*w_interp + Vsg*Vsg);
1917  const Real val_yz = std::sqrt(v_interp*v_interp + w_interp*w_interp + Vsg*Vsg);
1918  Gpu::deviceReduceSum(&plane_avg[iavg_xz], val_xz, handler);
1919  Gpu::deviceReduceSum(&plane_avg[iavg_yz], val_yz, handler);
1920  }
1921  });
1922  } else {
1923  auto k_arr = use_spatial_indices
1924  ? k_indx->const_array(mfi) : Array4<const int>{};
1925  auto j_arr = j_indx ? j_indx->const_array(mfi) : Array4<const int> {};
1926  auto i_arr = i_indx ? i_indx->const_array(mfi) : Array4<const int> {};
1927  const Box k_box = use_spatial_indices
1928  ? k_indx->fabbox(mfi.index()) : Box{};
1929  ParallelFor(Gpu::KernelInfo().setReduction(true), pbx, [=]
1930  AMREX_GPU_DEVICE(int i, int j, int k, Gpu::Handler const& handler) noexcept
1931  {
1932  const int ki = use_spatial_indices
1933  ? max(k_box.smallEnd(0), min(k_box.bigEnd(0), i)) : i;
1934  const int kj = use_spatial_indices
1935  ? max(k_box.smallEnd(1), min(k_box.bigEnd(1), j)) : j;
1936  const int kk = use_spatial_indices
1937  ? max(k_box.smallEnd(2), min(k_box.bigEnd(2), k)) : k;
1938  const int ref = use_spatial_indices
1939  ? k_arr(ki,kj,kk) : wall_normal_ref;
1940  int mi = i_arr ? i_arr(ki,kj,k) : i;
1941  int mj = j_arr ? j_arr(ki,kj,k) : j;
1942  int mk = k;
1943  if (dir == 0) {
1944  mi = ref;
1945  } else if (dir == 1) {
1946  mj = ref;
1947  } else {
1948  mk = ref;
1949  }
1950  const Real u_val = myhalf * (u_mf_arr(mi,mj,mk) + u_mf_arr(mi+1,mj ,mk));
1951  const Real v_val = myhalf * (v_mf_arr(mi,mj,mk) + v_mf_arr(mi ,mj+1,mk));
1952  const Real val = std::sqrt(u_val*u_val + v_val*v_val + Vsg*Vsg);
1953  Gpu::deviceReduceSum(&plane_avg[iavg], val, handler);
1954 
1955  if (dir < 2) {
1956  // averages for mean velocity on XZ and YZ planes
1957  const Real w_val = myhalf * (w_mf_arr(mi,mj,mk) + w_mf_arr(mi,mj,mk+1));
1958  const Real val_xz = std::sqrt(u_val*u_val + w_val*w_val + Vsg*Vsg);
1959  const Real val_yz = std::sqrt(v_val*v_val + w_val*w_val + Vsg*Vsg);
1960  Gpu::deviceReduceSum(&plane_avg[iavg_xz], val_xz, handler);
1961  Gpu::deviceReduceSum(&plane_avg[iavg_yz], val_yz, handler);
1962  }
1963  });
1964  }
1965  }
1966  }
1967 
1968  // Copy to host and sum across procs
1969  Gpu::copy(Gpu::deviceToHost, pavg.begin(), pavg.end(), plane_average.begin());
1970  ParallelDescriptor::ReduceRealSum(plane_average.data(), static_cast<int>(plane_average.size()));
1971 
1972  // No spatial variation with plane averages
1973  for (int iavg(0); iavg < m_navg; ++iavg){
1974  plane_average[iavg] *= denom[iavg]*d_fact_new;
1975  plane_average[iavg] += val_old[iavg];
1976  averages[iavg]->setVal(plane_average[iavg]);
1977  }
1978 }
dxInv
amrex::GpuArray< Real, AMREX_SPACEDIM > dxInv
Definition: ERF_InitCustomPertVels_ParticleTests.H:17
MOSTAverage::m_y_pos
amrex::Vector< std::unique_ptr< amrex::MultiFab > > m_y_pos
Definition: ERF_MOSTAverage.H:572
MOSTAverage::m_i_indx
amrex::Vector< std::unique_ptr< amrex::iMultiFab > > m_i_indx
Definition: ERF_MOSTAverage.H:574
MOSTAverage::m_z_phys_nd
amrex::Vector< amrex::MultiFab * > m_z_phys_nd
Definition: ERF_MOSTAverage.H:557
MOSTAverage::m_x_pos
amrex::Vector< std::unique_ptr< amrex::MultiFab > > m_x_pos
Definition: ERF_MOSTAverage.H:571
MOSTAverage::m_rot_fields
amrex::Vector< amrex::Vector< std::unique_ptr< amrex::MultiFab > > > m_rot_fields
Definition: ERF_MOSTAverage.H:579
MOSTAverage::m_z_pos
amrex::Vector< std::unique_ptr< amrex::MultiFab > > m_z_pos
Definition: ERF_MOSTAverage.H:573
MOSTAverage::m_ncell_plane
amrex::Vector< amrex::Vector< int > > m_ncell_plane
Definition: ERF_MOSTAverage.H:583
MOSTAverage::m_terrain_type
TerrainType m_terrain_type
Definition: ERF_MOSTAverage.H:560
MOSTAverage::m_norm_vec
bool m_norm_vec
Definition: ERF_MOSTAverage.H:595
MOSTAverage::m_j_indx
amrex::Vector< std::unique_ptr< amrex::iMultiFab > > m_j_indx
Definition: ERF_MOSTAverage.H:575
MOSTAverage::m_interp
bool m_interp
Definition: ERF_MOSTAverage.H:594
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:415
MOSTAverage::m_k_indx
amrex::Vector< std::unique_ptr< amrex::iMultiFab > > m_k_indx
Definition: ERF_MOSTAverage.H:576

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)

Fill averages for the point or region policy.

Parameters
[in]levlevel index

Function to compute average over local region.

Parameters
[in]levCurrent level
2088 {
2089  // Peel back the level
2090  auto& fields = m_fields[lev];
2091  auto& rot_fields = m_rot_fields[lev];
2092  auto& averages = m_averages[lev];
2093  const auto & geom = m_geom[lev];
2094 
2095  auto& z_phys = m_z_phys_nd[lev];
2096  auto& x_pos = m_x_pos[lev];
2097  auto& y_pos = m_y_pos[lev];
2098  auto& z_pos = m_z_pos[lev];
2099 
2100  auto& i_indx = m_i_indx[lev];
2101  auto& j_indx = m_j_indx[lev];
2102  auto& k_indx = m_k_indx[lev];
2103 
2104  const int dir = m_face.coordDir();
2105  // The non-terrain reference is constant on a level. Do not index the
2106  // cell-centered k-index FAB with an MFIter from a staggered velocity FAB;
2107  // distributed BoxArrays have the same box ordering but different valid
2108  // extents. Spatial index arrays are retained for terrain-normal cases.
2109  const bool fitted_terrain =
2110  (m_terrain_type == TerrainType::StaticFittedMesh) ||
2111  (m_terrain_type == TerrainType::MovingFittedMesh);
2112  const bool use_spatial_indices =
2113  m_norm_vec || (dir == 2 && fitted_terrain);
2114  if (!use_spatial_indices) {
2115  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
2116  m_k_indx[lev] != nullptr,
2117  "Region averaging requires a reference-index field.");
2118  }
2119  const int wall_normal_ref = use_spatial_indices ? 0 : m_k_indx[lev]->min(0);
2120 
2121  // Set factors for time averaging
2122  Real d_fact_new, d_fact_old;
2123  if (m_t_avg && m_t_init[lev]) {
2124  d_fact_new = m_fact_new;
2125  d_fact_old = m_fact_old;
2126  } else {
2127  d_fact_new = one;
2128  d_fact_old = zero;
2129  }
2130 
2131  // Number of cells contained in the local average
2132  const Real denom = one / (Real) m_ncell_region;
2133 
2134  // Capture radius for device
2135  int d_radius = m_radius;
2136 
2137  // NOTE: With more than one level we must explicitly fill the ghost cells of
2138  // the averages. An isolated fine patch has no neighboring box for the
2139  // FillBoundary below to communicate with, while the fields do carry
2140  // valid ghost data (filled from the coarse level), so we compute the
2141  // averages there as well. Where fine boxes do abut, the FillBoundary
2142  // overwrites what we compute with the neighbor's valid data.
2143  const IntVect ng_fill = get_ng_fill(lev);
2144 
2145  //
2146  //----------------------------------------------------------
2147  // Averages for U,V,W,T,Qv
2148  //----------------------------------------------------------
2149  //
2150  for (int imf(0); imf < 5; ++imf) {
2151 
2152  int sm_index = 0;
2153  if (m_face.isLow()) {
2154  sm_index = m_geom[lev].Domain().smallEnd(dir);
2155  } else {
2156  sm_index = m_geom[lev].Domain().bigEnd(dir);
2157  if (imf < 3 && imf == dir) {
2158  sm_index += 1;
2159  }
2160  }
2161 
2162  const int normal_face_offset =
2163  (!m_face.isLow() && imf < 3 && imf == dir) ? 1 : 0;
2164 
2165  // Continue if no valid Qv pointer
2166  if (!fields[imf]) continue;
2167 
2168 #ifdef _OPENMP
2169 #pragma omp parallel if (Gpu::notInLaunchRegion())
2170 #endif
2171  for (MFIter mfi(*fields[imf], TileNoZ()); mfi.isValid(); ++mfi) {
2172  Box pbx = mfi.growntilebox(ng_fill);
2173 
2174  if (m_face.isLow()) {
2175  if (mfi.validbox().smallEnd(dir) != sm_index ||
2176  pbx.smallEnd(dir) != sm_index) {
2177  continue;
2178  }
2179  } else {
2180  if (mfi.validbox().bigEnd(dir) != sm_index ||
2181  pbx.bigEnd(dir) != sm_index) {
2182  continue;
2183  }
2184  }
2185 
2186  pbx.setSmall(dir, sm_index); pbx.setBig(dir, sm_index);
2187 
2188  auto mf_arr = (m_rotate && imf != 2) ? rot_fields[imf]->const_array(mfi) :
2189  fields[imf]->const_array(mfi);
2190  auto ma_arr = averages[imf]->array(mfi);
2191 
2192  if (m_interp) {
2193  const auto plo = geom.ProbLoArray();
2194  const auto dx = geom.CellSizeArray();
2195  const auto dxInv = geom.InvCellSizeArray();
2196  const auto z_phys_arr = z_phys->const_array(mfi);
2197  auto x_pos_arr = x_pos->array(mfi);
2198  auto y_pos_arr = y_pos->array(mfi);
2199  auto z_pos_arr = z_pos->array(mfi);
2200  ParallelFor(pbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
2201  {
2202  ma_arr(i,j,k) *= d_fact_old;
2203 
2204  Real met_h_zeta = Compute_h_zeta_AtCellCenter(i,j,k,dxInv,z_phys_arr);
2205  for (int lk(-d_radius); lk <= (d_radius); ++lk) {
2206  for (int lj(-d_radius); lj <= (d_radius); ++lj) {
2207  for (int li(-d_radius); li <= (d_radius); ++li) {
2208  Real interp{0};
2209  Real xp = x_pos_arr(i+li,j+lj,k);
2210  Real yp = y_pos_arr(i+li,j+lj,k);
2211  Real zp = z_pos_arr(i+li,j+lj,k) + met_h_zeta*lk*dx[2];
2212  trilinear_interp_T(xp, yp, zp, &interp, mf_arr, z_phys_arr, plo, dxInv, 1);
2213  Real val = denom * interp * d_fact_new;
2214  ma_arr(i,j,k) += val;
2215  }
2216  }
2217  }
2218  });
2219  } else {
2220  auto k_arr = use_spatial_indices
2221  ? k_indx->const_array(mfi) : Array4<const int>{};
2222  auto j_arr = j_indx ? j_indx->const_array(mfi) : Array4<const int> {};
2223  auto i_arr = i_indx ? i_indx->const_array(mfi) : Array4<const int> {};
2224  const Box k_box = use_spatial_indices
2225  ? k_indx->fabbox(mfi.index()) : Box{};
2226  ParallelFor(pbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
2227  {
2228  const int ki = use_spatial_indices
2229  ? max(k_box.smallEnd(0), min(k_box.bigEnd(0), i)) : i;
2230  const int kj = use_spatial_indices
2231  ? max(k_box.smallEnd(1), min(k_box.bigEnd(1), j)) : j;
2232  const int kk = use_spatial_indices
2233  ? max(k_box.smallEnd(2), min(k_box.bigEnd(2), k)) : k;
2234  const int ref = (use_spatial_indices
2235  ? k_arr(ki,kj,kk) : wall_normal_ref) + normal_face_offset;
2236  int mi = i_arr ? i_arr(ki,kj,k) : i;
2237  int mj = j_arr ? j_arr(ki,kj,k) : j;
2238  int mk = k;
2239  if (dir == 0) {
2240  mi = ref;
2241  } else if (dir == 1) {
2242  mj = ref;
2243  } else {
2244  mk = ref;
2245  }
2246 
2247  ma_arr(i,j,k) *= d_fact_old;
2248  for (int lk(mk-d_radius); lk <= (mk+d_radius); ++lk) {
2249  for (int lj(mj-d_radius); lj <= (mj+d_radius); ++lj) {
2250  for (int li(mi-d_radius); li <= (mi+d_radius); ++li) {
2251  Real val = denom * mf_arr(li, lj, lk) * d_fact_new;
2252  ma_arr(i,j,k) += val;
2253  }
2254  }
2255  }
2256  });
2257  }
2258  } // MFiter
2259 
2260  // Fill the ghost cells we could not compute above
2261  //***********************************************************************************
2262  extrap_ghost_cells(lev,imf,ng_fill);
2263 
2264  // Fill interior ghost cells and any ghost cells outside a periodic domain
2265  //***********************************************************************************
2266  fill_planar_boundary(lev, *averages[imf]);
2267  } // imf
2268 
2269  //
2270  //----------------------------------------------------------
2271  // Averages for virtual potential temperature
2272  //----------------------------------------------------------
2273  //
2274  if (fields[4]) // We have water vapor
2275  {
2276  int sm_index = 0;
2277  if (m_face.isLow()) {
2278  sm_index = m_geom[lev].Domain().smallEnd(dir);
2279  } else {
2280  sm_index = m_geom[lev].Domain().bigEnd(dir);
2281  }
2282 
2283  int iavg = 5;
2284 
2285 #ifdef _OPENMP
2286 #pragma omp parallel if (Gpu::notInLaunchRegion())
2287 #endif
2288  for (MFIter mfi(*fields[4], TileNoZ()); mfi.isValid(); ++mfi) {
2289  Box pbx = mfi.growntilebox(ng_fill);
2290 
2291  if (m_face.isLow()) {
2292  if (mfi.validbox().smallEnd(dir) != sm_index ||
2293  pbx.smallEnd(dir) != sm_index) {
2294  continue;
2295  }
2296  } else {
2297  if (mfi.validbox().bigEnd(dir) != sm_index ||
2298  pbx.bigEnd(dir) != sm_index) {
2299  continue;
2300  }
2301  }
2302 
2303  pbx.setSmall(dir, sm_index); pbx.setBig(dir, sm_index);
2304 
2305  const Array4<Real const>& T_mf_arr = fields[3]->const_array(mfi);
2306  const Array4<Real const>& qv_mf_arr = (fields[4])? fields[4]->const_array(mfi) : Array4<const Real>{};
2307  const Array4<Real const>& qr_mf_arr = (fields[5])? fields[5]->const_array(mfi) : Array4<const Real>{};
2308  auto ma_arr = averages[iavg]->array(mfi);
2309 
2310  if (m_interp) {
2311  const auto plo = geom.ProbLoArray();
2312  const auto dx = geom.CellSizeArray();
2313  const auto dxInv = geom.InvCellSizeArray();
2314  const auto z_phys_arr = z_phys->const_array(mfi);
2315  auto x_pos_arr = x_pos->array(mfi);
2316  auto y_pos_arr = y_pos->array(mfi);
2317  auto z_pos_arr = z_pos->array(mfi);
2318  ParallelFor(pbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
2319  {
2320  ma_arr(i,j,k) *= d_fact_old;
2321 
2322  Real met_h_zeta = Compute_h_zeta_AtCellCenter(i,j,k,dxInv,z_phys_arr);
2323  for (int lk(-d_radius); lk <= (d_radius); ++lk) {
2324  for (int lj(-d_radius); lj <= (d_radius); ++lj) {
2325  for (int li(-d_radius); li <= (d_radius); ++li) {
2326  Real T_interp{0};
2327  Real qv_interp{0};
2328  Real xp = x_pos_arr(i+li,j+lj,k);
2329  Real yp = y_pos_arr(i+li,j+lj,k);
2330  Real zp = z_pos_arr(i+li,j+lj,k) + met_h_zeta*lk*dx[2];
2331  trilinear_interp_T(xp, yp, zp, &T_interp, T_mf_arr, z_phys_arr, plo, dxInv, 1);
2332  trilinear_interp_T(xp, yp, zp, &qv_interp, qv_mf_arr, z_phys_arr, plo, dxInv, 1);
2333  Real vfac;
2334  if (qr_mf_arr) {
2335  // We also have liquid water
2336  Real qr_interp{0};
2337  trilinear_interp_T(x_pos_arr(i,j,k), y_pos_arr(i,j,k), z_pos_arr(i,j,k),
2338  &qr_interp, qr_mf_arr, z_phys_arr, plo, dxInv, 1);
2339  vfac = one + epsv*qv_interp - qr_interp;
2340  } else {
2341  vfac = one + epsv*qv_interp;
2342  }
2343  const Real mag = T_interp * vfac;
2344  const Real val = denom * mag * d_fact_new;
2345  ma_arr(i,j,k) += val;
2346  }
2347  }
2348  }
2349  });
2350  } else {
2351  auto k_arr = use_spatial_indices
2352  ? k_indx->const_array(mfi) : Array4<const int>{};
2353  auto j_arr = j_indx ? j_indx->const_array(mfi) : Array4<const int> {};
2354  auto i_arr = i_indx ? i_indx->const_array(mfi) : Array4<const int> {};
2355  const Box k_box = use_spatial_indices
2356  ? k_indx->fabbox(mfi.index()) : Box{};
2357  ParallelFor(pbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
2358  {
2359  const int ki = use_spatial_indices
2360  ? max(k_box.smallEnd(0), min(k_box.bigEnd(0), i)) : i;
2361  const int kj = use_spatial_indices
2362  ? max(k_box.smallEnd(1), min(k_box.bigEnd(1), j)) : j;
2363  const int kk = use_spatial_indices
2364  ? max(k_box.smallEnd(2), min(k_box.bigEnd(2), k)) : k;
2365  const int ref = use_spatial_indices
2366  ? k_arr(ki,kj,kk) : wall_normal_ref;
2367  int mi = i_arr ? i_arr(ki,kj,k) : i;
2368  int mj = j_arr ? j_arr(ki,kj,k) : j;
2369  int mk = k;
2370  if (dir == 0) {
2371  mi = ref;
2372  } else if (dir == 1) {
2373  mj = ref;
2374  } else {
2375  mk = ref;
2376  }
2377 
2378  ma_arr(i,j,k) *= d_fact_old;
2379  for (int lk(mk-d_radius); lk <= (mk+d_radius); ++lk) {
2380  for (int lj(mj-d_radius); lj <= (mj+d_radius); ++lj) {
2381  for (int li(mi-d_radius); li <= (mi+d_radius); ++li) {
2382  Real vfac;
2383  if (qr_mf_arr) {
2384  // We also have liquid water
2385  vfac = one + epsv*qv_mf_arr(li,lj,lk) - qr_mf_arr(li,lj,lk);
2386  } else {
2387  vfac = one + epsv*qv_mf_arr(li,lj,lk);
2388  }
2389  const Real mag = T_mf_arr(li,lj,lk) * vfac;
2390  const Real val = denom * mag * d_fact_new;
2391  ma_arr(i,j,k) += val;
2392  }
2393  }
2394  }
2395  });
2396  }
2397  } // MFiter
2398 
2399  // Fill the ghost cells we could not compute above
2400  //***********************************************************************************
2401  extrap_ghost_cells(lev,iavg,ng_fill);
2402 
2403  // Fill interior ghost cells and any ghost cells outside a periodic domain
2404  //***********************************************************************************
2405  fill_planar_boundary(lev, *averages[iavg]);
2406 
2407  }
2408  else // copy temperature
2409  {
2410  int iavg = m_navg - 4;
2411  IntVect ng = averages[iavg]->nGrowVect();
2412  MultiFab::Copy(*(averages[iavg]),*(averages[3]),0,0,1,ng);
2413  }
2414 
2415  //
2416  //----------------------------------------------------------
2417  // Averages for the tangential velocity magnitude
2418  //----------------------------------------------------------
2419  //
2420  {
2421  int sm_index = 0;
2422  if (m_face.isLow()) {
2423  sm_index = m_geom[lev].Domain().smallEnd(dir);
2424  } else {
2425  sm_index = m_geom[lev].Domain().bigEnd(dir);
2426  }
2427 
2428  const int imf_cc = 3;
2429  const int imf = 0;
2430  const int iavg = m_navg - 3;
2431  const int iavg_xz = m_navg - 2;
2432  const int iavg_yz = m_navg - 1;
2433 
2434  const Real Vsg = m_Vsg[lev];
2435 
2436 #ifdef _OPENMP
2437 #pragma omp parallel if (Gpu::notInLaunchRegion())
2438 #endif
2439  for (MFIter mfi(*fields[imf_cc], TileNoZ()); mfi.isValid(); ++mfi) {
2440  Box pbx = mfi.growntilebox(ng_fill);
2441 
2442  if (m_face.isLow()) {
2443  if (mfi.validbox().smallEnd(dir) != sm_index ||
2444  pbx.smallEnd(dir) != sm_index) {
2445  continue;
2446  }
2447  } else {
2448  if (mfi.validbox().bigEnd(dir) != sm_index ||
2449  pbx.bigEnd(dir) != sm_index) {
2450  continue;
2451  }
2452  }
2453 
2454  pbx.setSmall(dir, sm_index); pbx.setBig(dir, sm_index);
2455 
2456  auto u_mf_arr = (m_rotate) ? rot_fields[imf ]->const_array(mfi) :
2457  fields[imf ]->const_array(mfi);
2458  auto v_mf_arr = (m_rotate) ? rot_fields[imf+1]->const_array(mfi) :
2459  fields[imf+1]->const_array(mfi);
2460  auto w_mf_arr = fields[imf+2]->const_array(mfi);
2461  auto ma_arr = averages[iavg]->array(mfi);
2462  auto ma_xz_arr = averages[iavg_xz]->array(mfi);
2463  auto ma_yz_arr = averages[iavg_yz]->array(mfi);
2464 
2465 
2466  if (m_interp) {
2467  const auto plo = geom.ProbLoArray();
2468  const auto dx = geom.CellSizeArray();
2469  const auto dxInv = geom.InvCellSizeArray();
2470  const auto z_phys_arr = z_phys->const_array(mfi);
2471  auto x_pos_arr = x_pos->array(mfi);
2472  auto y_pos_arr = y_pos->array(mfi);
2473  auto z_pos_arr = z_pos->array(mfi);
2474  ParallelFor(pbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
2475  {
2476  ma_arr(i,j,k) *= d_fact_old;
2477 
2478  Real met_h_zeta = Compute_h_zeta_AtCellCenter(i,j,k,dxInv,z_phys_arr);
2479  for (int lk(-d_radius); lk <= (d_radius); ++lk) {
2480  for (int lj(-d_radius); lj <= (d_radius); ++lj) {
2481  for (int li(-d_radius); li <= (d_radius); ++li) {
2482  Real u_interp{0};
2483  Real v_interp{0};
2484  Real xp = x_pos_arr(i+li,j+lj,k);
2485  Real yp = y_pos_arr(i+li,j+lj,k);
2486  Real zp = z_pos_arr(i+li,j+lj,k) + met_h_zeta*lk*dx[2];
2487  trilinear_interp_T(xp, yp, zp, &u_interp, u_mf_arr, z_phys_arr, plo, dxInv, 1);
2488  trilinear_interp_T(xp, yp, zp, &v_interp, v_mf_arr, z_phys_arr, plo, dxInv, 1);
2489  const Real mag = std::sqrt(u_interp*u_interp + v_interp*v_interp + Vsg*Vsg);
2490  Real val = denom * mag * d_fact_new;
2491  ma_arr(i,j,k) += val;
2492  }
2493  }
2494  }
2495  });
2496  } else {
2497  auto k_arr = use_spatial_indices
2498  ? k_indx->const_array(mfi) : Array4<const int>{};
2499  auto j_arr = j_indx ? j_indx->const_array(mfi) : Array4<const int> {};
2500  auto i_arr = i_indx ? i_indx->const_array(mfi) : Array4<const int> {};
2501  const Box k_box = use_spatial_indices
2502  ? k_indx->fabbox(mfi.index()) : Box{};
2503  ParallelFor(pbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
2504  {
2505  const int ki = use_spatial_indices
2506  ? max(k_box.smallEnd(0), min(k_box.bigEnd(0), i)) : i;
2507  const int kj = use_spatial_indices
2508  ? max(k_box.smallEnd(1), min(k_box.bigEnd(1), j)) : j;
2509  const int kk = use_spatial_indices
2510  ? max(k_box.smallEnd(2), min(k_box.bigEnd(2), k)) : k;
2511  const int ref = use_spatial_indices
2512  ? k_arr(ki,kj,kk) : wall_normal_ref;
2513  int mi = i_arr ? i_arr(ki,kj,k) : i;
2514  int mj = j_arr ? j_arr(ki,kj,k) : j;
2515  int mk = k;
2516  if (dir == 0) {
2517  mi = ref;
2518  } else if (dir == 1) {
2519  mj = ref;
2520  } else {
2521  mk = ref;
2522  }
2523 
2524  ma_arr(i,j,k) *= d_fact_old;
2525 
2526  if (dir < 2) {
2527  ma_xz_arr(i,j,k) *= d_fact_old;
2528  ma_yz_arr(i,j,k) *= d_fact_old;
2529  } else {
2530  ma_xz_arr(i,j,k) = 0.0;
2531  ma_yz_arr(i,j,k) = 0.0;
2532  }
2533  for (int lk(mk-d_radius); lk <= (mk+d_radius); ++lk) {
2534  for (int lj(mj-d_radius); lj <= (mj+d_radius); ++lj) {
2535  for (int li(mi-d_radius); li <= (mi+d_radius); ++li) {
2536  const Real u_val = myhalf * (u_mf_arr(li,lj,lk) + u_mf_arr(li+1,lj ,lk));
2537  const Real v_val = myhalf * (v_mf_arr(li,lj,lk) + v_mf_arr(li ,lj+1,lk));
2538  const Real mag = std::sqrt(u_val*u_val + v_val*v_val + Vsg*Vsg);
2539  Real val = denom * mag * d_fact_new;
2540  ma_arr(i,j,k) += val;
2541 
2542  if (dir < 2) {
2543  // averages for mean velocity on XZ and YZ planes
2544  const Real w_val = myhalf * (w_mf_arr(li,lj,lk) + w_mf_arr(li,lj,lk+1));
2545  const Real val_xz = std::sqrt(u_val*u_val + w_val*w_val + Vsg*Vsg);
2546  const Real val_yz = std::sqrt(v_val*v_val + w_val*w_val + Vsg*Vsg);
2547 
2548  ma_xz_arr(i,j,k) += denom * val_xz * d_fact_new;
2549  ma_yz_arr(i,j,k) += denom * val_yz * d_fact_new;
2550  }
2551  }
2552  }
2553  }
2554  });
2555  }
2556  } // MFiter
2557 
2558  // Fill the ghost cells we could not compute above
2559  //***********************************************************************************
2560  extrap_ghost_cells(lev,iavg,ng_fill);
2561 
2562  // Fill interior ghost cells and any ghost cells outside a periodic domain
2563  //***********************************************************************************
2564  fill_planar_boundary(lev, *averages[iavg]);
2565  if (dir < 2) {
2566  fill_planar_boundary(lev, *averages[iavg_xz]);
2567  fill_planar_boundary(lev, *averages[iavg_yz]);
2568  }
2569  }
2570 
2571  // NOTE: Checking periodicity with the geom structure is not
2572  // sufficient at higher levels. The BA may be contained
2573  // within the domain and it's exterior ghost cells filled
2574  // from interpolation; yet the domain BCs are periodic.
2575 
2576  // Need to fill ghost cells outside the domain if not periodic
2577  bool not_per_x = !(geom.periodicity().isPeriodic(0));
2578  bool not_per_y = !(geom.periodicity().isPeriodic(1));
2579  const bool per_x = geom.periodicity().isPeriodic(0);
2580  const bool per_y = geom.periodicity().isPeriodic(1);
2581  const bool per_z = geom.periodicity().isPeriodic(2);
2582  Box cc_bnd_bx = (m_fields[lev][3]->boxArray()).minimalBox();
2583  Box domain = geom.Domain();
2584 
2585  if (domain.contains(cc_bnd_bx) || (not_per_x || not_per_y)) {
2586  for (int iavg(0); iavg < m_navg; ++iavg) {
2587  IntVect ng = averages[iavg]->nGrowVect();
2588 
2589  // NOTE: Level 0 spans the whole domain, but finer
2590  // levels do not have such a restriction.
2591  // For now, use the bounding box of the boxArray.
2592 
2593  // NOTE: The fields and averages have different indexing.
2594  // The averages are: U/V/W/T/Qv/Tv/Umag_XY/Umag_XZ/Umag_YZ
2595  // The fields are: U/V/W/T/Qv/Qr
2596  // We clip iavg at 3 since all the remaining data is CC
2597 
2598  int imf = min(iavg,3);
2599 
2600  // For a z face, tangential extrapolation follows the level's
2601  // actual BoxArray extent. Fine patches may occupy only part of
2602  // the domain, and their interior ghosts must remain untouched.
2603  Box bnd_bx = m_fields[lev][imf]->boxArray().minimalBox();
2604 
2605  // Domain box in the index type of the current source field.
2606  // This is the physical-domain skip extent, including the
2607  // high-side nodal offset for velocity fields.
2608  Box dom_bx = convert(domain, fields[imf]->boxArray().ixType());
2609 
2610  int sm_index = 0;
2611  if (m_face.isLow()) {
2612  sm_index = dom_bx.smallEnd(dir);
2613  } else {
2614  sm_index = dom_bx.bigEnd(dir);
2615  }
2616 #ifdef _OPENMP
2617 #pragma omp parallel if (Gpu::notInLaunchRegion())
2618 #endif
2619  for (MFIter mfi(*fields[imf], TileNoZ()); mfi.isValid(); ++mfi) {
2620  Box vbx = mfi.validbox();
2621  Box pbx = mfi.tilebox();
2622  Box gpbx = mfi.growntilebox(ng);
2623 
2624  if (dom_bx.contains(gpbx)) {
2625  continue;
2626  }
2627 
2628  if (m_face.isLow()) {
2629  if (vbx.smallEnd(dir) != sm_index ||
2630  pbx.smallEnd(dir) != sm_index) {
2631  continue;
2632  }
2633  gpbx.setBig(dir, sm_index);
2634  } else {
2635  if (vbx.bigEnd(dir) != sm_index ||
2636  pbx.bigEnd(dir) != sm_index) {
2637  continue;
2638  }
2639  gpbx.setSmall(dir, sm_index);
2640  }
2641 
2642  auto ma_arr = averages[iavg]->array(mfi);
2643 
2644  if (dir == 0) {
2645  int j_lo = dom_bx.smallEnd(1); int j_hi = dom_bx.bigEnd(1);
2646  int k_lo = dom_bx.smallEnd(2); int k_hi = dom_bx.bigEnd(2);
2647 
2648  ParallelFor(gpbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
2649  {
2650  // Periodic tangential ghosts were populated by the
2651  // FillBoundary above. Leave them untouched; this
2652  // kernel only supplies non-periodic ghost values.
2653  if ((per_y && (j < j_lo || j > j_hi)) ||
2654  (per_z && (k < k_lo || k > k_hi))) {
2655  return;
2656  }
2657 
2658  int lj, lk;
2659  lj = j < j_lo ? j_lo : j;
2660  lj = lj > j_hi ? j_hi : lj;
2661  lk = k < k_lo ? k_lo : k;
2662  lk = lk > k_hi ? k_hi : lk;
2663 
2664  ma_arr(i,j,k) = ma_arr(sm_index,lj,lk);
2665  });
2666  } else if (dir == 1) {
2667  int i_lo = dom_bx.smallEnd(0); int i_hi = dom_bx.bigEnd(0);
2668  int k_lo = dom_bx.smallEnd(2); int k_hi = dom_bx.bigEnd(2);
2669  ParallelFor(gpbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
2670  {
2671  // Periodic tangential ghosts were populated by the
2672  // FillBoundary above. Leave them untouched; this
2673  // kernel only supplies non-periodic ghost values.
2674  if ((per_x && (i < i_lo || i > i_hi)) ||
2675  (per_z && (k < k_lo || k > k_hi))) {
2676  return;
2677  }
2678 
2679  int li, lk;
2680  li = i < i_lo ? i_lo : i;
2681  li = li > i_hi ? i_hi : li;
2682  lk = k < k_lo ? k_lo : k;
2683  lk = lk > k_hi ? k_hi : lk;
2684 
2685  ma_arr(i,j,k) = ma_arr(li,sm_index,lk);
2686  });
2687  } else {
2688  int i_lo = bnd_bx.smallEnd(0); int i_hi = bnd_bx.bigEnd(0);
2689  int j_lo = bnd_bx.smallEnd(1); int j_hi = bnd_bx.bigEnd(1);
2690  ParallelFor(gpbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
2691  {
2692  // Periodic tangential ghosts were populated by the
2693  // FillBoundary above. Leave them untouched; this
2694  // kernel only supplies non-periodic ghost values.
2695  if ((per_x && (i < i_lo || i > i_hi)) ||
2696  (per_y && (j < j_lo || j > j_hi))) {
2697  return;
2698  }
2699 
2700  int li, lj;
2701  li = i < i_lo ? i_lo : i;
2702  li = li > i_hi ? i_hi : li;
2703  lj = j < j_lo ? j_lo : j;
2704  lj = lj > j_hi ? j_hi : lj;
2705 
2706  ma_arr(i,j,k) = ma_arr(li,lj,sm_index);
2707  });
2708  }
2709  } // MFiter
2710  } // iavg
2711  } // Not periodic
2712 }
Compute_h_zeta_AtCellCenter
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE 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:201
MOSTAverage::m_radius
int m_radius
Definition: ERF_MOSTAverage.H:588
MOSTAverage::fill_planar_boundary
void fill_planar_boundary(const int &lev, amrex::MultiFab &mf)
Definition: ERF_MOSTAverage.cpp:2069
MOSTAverage::m_ncell_region
int m_ncell_region
Definition: ERF_MOSTAverage.H:589
MOSTAverage::extrap_ghost_cells
void extrap_ghost_cells(const int &lev, const int &iavg, const amrex::IntVect &ng_fill)
Definition: ERF_MOSTAverage.cpp:1996
MOSTAverage::get_ng_fill
amrex::IntVect get_ng_fill(const int &lev) const
Definition: ERF_MOSTAverage.cpp:473
MicVar_Morr::ng
@ ng
Definition: ERF_Morrison.H:51

Referenced by compute_averages().

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

◆ do_time_averaging()

bool MOSTAverage::do_time_averaging ( ) const
inline

Return whether the averages are filtered in time (erf.most.time_average).

265 { return m_t_avg; }

Referenced by SurfaceLayer::mac_avg_is_time_averaged().

Here is the caller graph for this function:

◆ extrap_ghost_cells()

void MOSTAverage::extrap_ghost_cells ( const int &  lev,
const int &  iavg,
const amrex::IntVect &  ng_fill 
)

Fill the ghost cells of one average that cannot be computed from the field data with a zeroth-order extrapolation of the nearest computed value.

The averages hold as many ghost cells as the fields, which is more than can be computed from them (see get_ng_fill), so the outermost layers are filled here to leave no average undefined. This is done before the FillBoundary, which then overwrites whatever is shared with a neighbor.

Parameters
[in]levlevel index
[in]iavgaverage component index
[in]ng_fillnumber of ghost cells that hold computed data

Function to fill the ghost cells of one average that cannot be computed from the field data with a zeroth-order extrapolation of the nearest computed value.

The averages hold as many ghost cells as the fields, which is more than can be computed from them (see get_ng_fill), so the outermost layers are filled here to leave no average undefined. This is done before the FillBoundary, which then overwrites whatever is shared with a neighboring box.

Parameters
[in]levCurrent level
[in]iavgAverage component
[in]ng_fillNumber of ghost cells holding computed data
1999 {
2000  // Peel back the level
2001  auto& averages = m_averages[lev];
2002 
2003  // Extrapolation is performed directly in each average's own index space;
2004  // the source fields use different centering and are not needed here.
2005 
2006  IntVect ng = averages[iavg]->nGrowVect();
2007  // Only the selected wall-normal direction is collapsed. For x/y faces
2008  // the z ghosts are tangential and must be extrapolated/fillable.
2009  ng[m_face.coordDir()] = 0;
2010 
2011  // Everything we hold was computed above
2012  if (ng.allLE(ng_fill)) { return; }
2013 
2014  IntVect ng_computed = ng_fill;
2015  ng_computed[m_face.coordDir()] = 0;
2016 
2017 #ifdef _OPENMP
2018 #pragma omp parallel if (Gpu::notInLaunchRegion())
2019 #endif
2020  for (MFIter mfi(*averages[iavg], false); mfi.isValid(); ++mfi) {
2021  // Use the average's index space. The source field and the average
2022  // have different layouts for U/V/W and for the selected wall face.
2023  const Box avg_fab = averages[iavg]->fabbox(mfi.index());
2024  Box gpbx = averages[iavg]->boxArray()[mfi.index()];
2025  gpbx.grow(ng);
2026  gpbx &= avg_fab;
2027 
2028  // Region of this average FAB that was computed from field data.
2029  Box cbx = averages[iavg]->boxArray()[mfi.index()];
2030  cbx.grow(ng_computed);
2031  cbx &= avg_fab;
2032 
2033  if (cbx.contains(gpbx)) { continue; }
2034 
2035  auto ma_arr = averages[iavg]->array(mfi);
2036 
2037  int i_lo = cbx.smallEnd(0); int i_hi = cbx.bigEnd(0);
2038  int j_lo = cbx.smallEnd(1); int j_hi = cbx.bigEnd(1);
2039  int k_lo = cbx.smallEnd(2); int k_hi = cbx.bigEnd(2);
2040  ParallelFor(gpbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
2041  {
2042  int li, lj, lk;
2043  li = i < i_lo ? i_lo : i;
2044  li = li > i_hi ? i_hi : li;
2045  lj = j < j_lo ? j_lo : j;
2046  lj = lj > j_hi ? j_hi : lj;
2047  lk = k < k_lo ? k_lo : k;
2048  lk = lk > k_hi ? k_hi : lk;
2049 
2050  ma_arr(i,j,k) = ma_arr(li,lj,lk);
2051  });
2052  } // MFiter
2053 }

Referenced by compute_region_averages().

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

◆ fill_planar_boundary()

void MOSTAverage::fill_planar_boundary ( const int &  lev,
amrex::MultiFab &  mf 
)

Fill the ghost cells of one planar average, and the valid region of its uncomputed copies when the 3D BoxArray is split in z (see PlanarBoundary).

Parameters
[in]levlevel index
[in,out]mfplanar average to fill

Function to fill the ghost cells of one planar average.

The averages hold one box per 3D box, so a 3D BoxArray split in the face-normal direction gives duplicate planar boxes, of which only the surface copy is computed by compute_region_averages; a FillBoundary could then fill a ghost cell from an uncomputed copy (see PlanarBoundary). With the split, the valid region of the uncomputed copies is filled as well. With EB terrain the averages are computed on every box and FillBoundary is well defined.

Parameters
[in]levCurrent level
[in,out]mfPlanar average to fill
2070 {
2071  const int dir = m_face.coordDir();
2072  const Periodicity period = tangential_periodicity(m_geom[lev], dir);
2073 
2074  if (m_terrain_type == TerrainType::EB) {
2075  mf.FillBoundary(mf.nGrowVect(), period);
2076  } else {
2077  m_planar_bndry[lev].fill(mf, period);
2078  }
2079 }
MOSTAverage::m_planar_bndry
amrex::Vector< PlanarBoundary > m_planar_bndry
Definition: ERF_MOSTAverage.H:578

Referenced by compute_region_averages().

Here is the caller graph for this function:

◆ get_average() [1/2]

amrex::MultiFab* MOSTAverage::get_average ( const int &  lev,
const int &  comp 
)
inline

Return one 2D average MultiFab for modification (used to restore the exponential filter state from a checkpoint).

Parameters
[in]levlevel index
[in]compaverage component index
260 { return m_averages[lev][comp].get(); }

◆ get_average() [2/2]

const amrex::MultiFab* MOSTAverage::get_average ( const int &  lev,
const int &  comp 
) const
inline

Return one 2D average MultiFab.

Parameters
[in]levlevel index
[in]compaverage component index
251 { return m_averages[lev][comp].get(); }

Referenced by SurfaceLayer::get_mac_avg(), and SurfaceLayer::get_mac_avg_ptr().

Here is the caller graph for this function:

◆ get_k_indices()

const amrex::iMultiFab* MOSTAverage::get_k_indices ( const int &  lev) const
inline

Return the k-index iMultiFab.

Parameters
[in]levlevel index
336 { return m_k_indx[lev].get(); }

◆ get_navg()

int MOSTAverage::get_navg ( ) const
inline

Return the number of average components.

270 { return m_navg; }

Referenced by SurfaceLayer::get_num_mac_avg().

Here is the caller graph for this function:

◆ get_ng_fill()

IntVect MOSTAverage::get_ng_fill ( const int &  lev) const

Return the number of ghost cells of the 2D data (averages, indices and positions) that can be filled from the field data at this level.

With more than one level we must explicitly fill the ghost cells of the 2D data: an isolated fine patch has no neighboring box for a FillBoundary to communicate with, while the fields themselves do carry valid ghost data (filled from the coarse level). The 2D data are allocated with as many ghost cells as the fields they are built from, but an average reads the fields over m_radius cells in each direction and the tangential velocity magnitude (like the interpolation stencil) reaches one cell beyond that. So a ghost cell of the 2D data can only be computed where the fields carry m_radius+1 ghost cells past it.

Parameters
[in]levlevel index

Function to return the number of ghost cells of the 2D data (averages, indices and positions) that can be filled from the field data.

With more than one level we must explicitly fill the ghost cells of the 2D data: an isolated fine patch has no neighboring box for a FillBoundary to communicate with, while the fields themselves do carry valid ghost data (filled from the coarse level). The 2D data are allocated with as many ghost cells as the fields they are built from, but an average reads the fields over m_radius cells in each direction and the tangential velocity magnitude (like the interpolation stencil) reaches one cell beyond that. So a ghost cell of the 2D data can only be computed where the fields carry m_radius+1 ghost cells past it.

Parameters
[in]levCurrent level
474 {
475  // The averages can hold no more than they were allocated with. Keep a
476  // separate limit for each direction: z is a tangential direction for
477  // x/y faces, while x/y are tangential for a z face.
478  IntVect ng_min = m_averages[lev][0]->nGrowVect();
479  for (int iavg(0); iavg < m_navg; ++iavg) {
480  const IntVect ng = m_averages[lev][iavg]->nGrowVect();
481  for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
482  ng_min[idim] = min(ng_min[idim], ng[idim]);
483  }
484  }
485 
486  // The velocities carry one fewer ghost cell than the CC fields, so they
487  // are what limits us in practice. Apply the same per-direction limit.
488  for (int imf(0); imf < m_nvar; ++imf) {
489  if (!m_fields[lev][imf]) { continue; }
490  const IntVect ng = m_fields[lev][imf]->nGrowVect();
491  for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
492  ng_min[idim] = min(ng_min[idim], ng[idim]);
493  }
494  }
495 
496  IntVect ng_fill(0);
497  for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
498  if (idim == m_face.coordDir()) { continue; }
499  // A radius-r average reads r cells in either direction. Staggered
500  // velocity interpolation reaches one additional cell on the high
501  // side, so that cell must also be available in the source field.
502  ng_fill[idim] = max(0, ng_min[idim] - (m_radius + 1));
503  }
504  return ng_fill;
505 }
MOSTAverage::m_nvar
int m_nvar
Definition: ERF_MOSTAverage.H:565

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

Here is the caller graph for this function:

◆ get_plane_average()

amrex::Vector<amrex::Real> MOSTAverage::get_plane_average ( const int &  lev) const
inline

Return the filtered plane averages, which hold the filter state for the plane and EB policies. The region policy holds its state in m_averages instead and never allocates these, in which case this is empty.

Parameters
[in]levlevel index
302  {
303  if (lev < static_cast<int>(m_plane_average.size())) { return m_plane_average[lev]; }
304  return amrex::Vector<amrex::Real>{};
305  }

Referenced by SurfaceLayer::get_mac_plane_avg().

Here is the caller graph for this function:

◆ get_zref()

amrex::MultiFab* MOSTAverage::get_zref ( const int &  lev) const
inline

Return z_ref, which may be computed from a specified k-index.

Parameters
[in]levlevel index
329 { return m_zref[lev].get(); }
MOSTAverage::m_zref
amrex::Vector< std::unique_ptr< amrex::MultiFab > > m_zref
Definition: ERF_MOSTAverage.H:570

Referenced by SurfaceLayer::get_zref(), and SurfaceLayer::get_zref_min().

Here is the caller graph for this function:

◆ in_cell_z()

AMREX_GPU_HOST_DEVICE static AMREX_INLINE bool MOSTAverage::in_cell_z ( const amrex::Real &  z,
const amrex::Real &  z_lo,
const amrex::Real &  z_hi 
)
inlinestatic

Whether a height lies in the vertical extent [z_lo, z_hi) of a cell.

The interval is closed below so that a query height landing exactly on a face (e.g. the default 10 m on a flat mesh with a 10 m first cell) belongs to the cell above it rather than to neither cell.

Parameters
[in]zquery height
[in]z_loheight of the lower face
[in]z_hiheight of the upper face
353  {
354  return ( (z >= z_lo) && (z < z_hi) );
355  }
Coord::z
@ z

Referenced by set_norm_indices_T(), and trilinear_interp_T().

Here is the caller graph for this function:

◆ k_index_below()

static int MOSTAverage::k_index_below ( const amrex::Vector< amrex::Real > &  zlevels,
const amrex::Real &  zref 
)
inlinestatic

Index of the highest cell whose center lies at or below a height above the bottom of the domain, on a column of staggered (face) heights. This is the stretched-mesh counterpart of floor(zref/dz - 1/2) on a uniform mesh.

Parameters
[in]zlevelsstaggered heights of the column (nz+1 faces)
[in]zrefheight above zlevels[0]
Returns
the cell index, or -1 if zref lies below the first cell center
368  {
369  const int nz = static_cast<int>(zlevels.size()) - 1;
370  int lk = -1;
371  for (int k(0); k < nz; ++k) {
372  if (cell_center_height(zlevels, k) <= zref) { lk = k; } else { break; }
373  }
374  return lk;
375  }
MOSTAverage::cell_center_height
static amrex::Real cell_center_height(const amrex::Vector< amrex::Real > &zlevels, const int &k)
Definition: ERF_MOSTAverage.H:384

Referenced by set_k_indices_N().

Here is the call graph for this function:
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 
)

Make MOST-average data structures at one level.

Parameters
[in]levlevel index
[in]vars_oldold-time state and velocity fields
[in]Theta_primprimitive potential-temperature field
[in]Qv_primprimitive water-vapor field
[in]Qr_primprimitive rain-water field
[in]z_phys_ndnodal physical-height field

Make MOSTAverage data structures at one level.

Parameters
[in]levCurrent level
[in]vars_oldState data used by the average calculator
[in]Theta_primPrimitive theta component at this level
[in]Qv_primPrimitive water-vapor component at this level
[in]Qr_primPrimitive rain-water component at this level
[in]z_phys_ndNodal physical height at this level
178 {
179  m_fields[lev].resize(m_nvar);
180  m_rot_fields[lev].resize(m_nvar-1);
181  m_averages[lev].resize(m_navg);
182  m_z_phys_nd[lev] = z_phys_nd.get();
183 
184  bool use_terrain_fitted_coords = ( (m_terrain_type == TerrainType::StaticFittedMesh) ||
185  (m_terrain_type == TerrainType::MovingFittedMesh) );
186 
187  const int dir = m_face.coordDir();
188 
189  bool use_eb = (m_terrain_type == TerrainType::EB);
190 
191  { // Nodal in x
192  auto& mf = *vars_old[Vars::xvel];
193  // Create a 2D ba, dm, & ghost cells
194  int sm_index;
195  const BoxArray& ba = mf.boxArray();
196  const Box dom_bx = convert(m_geom[lev].Domain(), ba.ixType());
197  if (m_face.isLow()) {
198  sm_index = dom_bx.smallEnd(dir);
199  } else {
200  sm_index = dom_bx.bigEnd(dir);
201  }
202 
203  BoxList bl2d = ba.boxList();
204  for (auto& b : bl2d) { b.setRange(dir, sm_index); }
205  BoxArray ba2d(std::move(bl2d));
206  const DistributionMapping& dm = mf.DistributionMap();
207  const int ncomp = 1;
208  IntVect ng = mf.nGrowVect();
209  ng[dir] = 0;
210 
211  m_fields[lev][0] = vars_old[Vars::xvel];
212  m_averages[lev][0] = std::make_unique<MultiFab>(ba2d,dm,ncomp,ng);
213  m_averages[lev][0]->setVal(bogus_large_value);
214  if (m_rotate) {
215  m_rot_fields[lev][0] = std::make_unique<MultiFab>(ba,dm,ncomp,ng);
216  MultiFab::Copy(*m_rot_fields[lev][0],mf,0,0,1,ng);
217  } else {
218  m_rot_fields[lev][0] = nullptr;
219  }
220  }
221  { // Nodal in y
222  auto& mf = *vars_old[Vars::yvel];
223  // Create a 2D ba, dm, & ghost cells
224  int sm_index;
225  const BoxArray& ba = mf.boxArray();
226  const Box dom_bx = convert(m_geom[lev].Domain(), ba.ixType());
227 
228  if (m_face.isLow()) {
229  sm_index = dom_bx.smallEnd(dir);
230  } else {
231  sm_index = dom_bx.bigEnd(dir);
232  }
233 
234  BoxList bl2d = ba.boxList();
235  for (auto& b : bl2d) { b.setRange(dir, sm_index); }
236  BoxArray ba2d(std::move(bl2d));
237  const DistributionMapping& dm = mf.DistributionMap();
238  const int ncomp = 1;
239  IntVect ng = mf.nGrowVect();
240  ng[dir] = 0;
241 
242  m_fields[lev][1] = vars_old[Vars::yvel];
243  m_averages[lev][1] = std::make_unique<MultiFab>(ba2d,dm,ncomp,ng);
244  m_averages[lev][1]->setVal(bogus_large_value);
245  if (m_rotate) {
246  m_rot_fields[lev][1] = std::make_unique<MultiFab>(ba,dm,ncomp,ng);
247  MultiFab::Copy(*m_rot_fields[lev][1],mf,0,0,1,ng);
248  } else {
249  m_rot_fields[lev][1] = nullptr;
250  }
251  }
252  { // Nodal in z
253  auto& mf = *vars_old[Vars::zvel];
254  // Create a 2D ba, dm, & ghost cells
255  int sm_index;
256  const BoxArray& ba = mf.boxArray();
257  const Box dom_bx = convert(m_geom[lev].Domain(), ba.ixType());
258 
259  if (m_face.isLow()) {
260  sm_index = dom_bx.smallEnd(dir);
261  } else {
262  sm_index = dom_bx.bigEnd(dir);
263  }
264 
265  BoxList bl2d = ba.boxList();
266  for (auto& b : bl2d) b.setRange(dir, sm_index);
267  BoxArray ba2d(std::move(bl2d));
268  const DistributionMapping& dm = mf.DistributionMap();
269  const int ncomp = 1;
270  IntVect ng = mf.nGrowVect();
271  ng[dir] = 0;
272 
273  m_fields[lev][2] = vars_old[Vars::zvel];
274  m_averages[lev][2] = std::make_unique<MultiFab>(ba2d,dm,ncomp,ng);
275  m_averages[lev][2]->setVal(1.E34);
276  }
277  { // CC vars
278  auto& mf = *Theta_prim;
279  // Create a 2D ba, dm, & ghost cells
280  int sm_index;
281  if (m_face.isLow()) {
282  sm_index = m_geom[lev].Domain().smallEnd(dir);
283  } else {
284  sm_index = m_geom[lev].Domain().bigEnd(dir);
285  }
286 
287  const BoxArray& ba = mf.boxArray();
288  BoxList bl2d = ba.boxList();
289  for (auto& b : bl2d) { b.setRange(dir, sm_index); }
290  BoxArray ba2d(std::move(bl2d));
291  const DistributionMapping& dm = mf.DistributionMap();
292  const int ncomp = 1;
293  const int incomp = 1;
294  IntVect ng = mf.nGrowVect();
295  ng[dir] = 0;
296 
297  // Get field pointers
298  m_fields[lev][3] = Theta_prim.get();
299  m_fields[lev][4] = Qv_prim.get();
300  m_fields[lev][5] = Qr_prim.get();
301 
302  if (!use_eb) {
303  const int surface_index = m_face.isLow()
304  ? m_geom[lev].Domain().smallEnd(dir)
305  : m_geom[lev].Domain().bigEnd(dir);
306  m_planar_bndry[lev].define(ba, ba2d, dm, surface_index, m_face.isLow(), dir);
307  }
308 
309  // Initialize remaining multifabs
310  for (int iavg(3); iavg < m_navg; ++iavg) {
311  m_averages[lev][iavg] = std::make_unique<MultiFab>(ba2d,dm,ncomp,ng);
312  m_averages[lev][iavg]->setVal(bogus_large_value);
313  }
314 
315  // Default to dry
316  m_averages[lev][4]->setVal(0.0); // Qv
317 
318  if (m_rotate) {
319  m_rot_fields[lev][3] = std::make_unique<MultiFab>(ba,dm,ncomp,ng);
320  m_rot_fields[lev][4] = std::make_unique<MultiFab>(ba,dm,ncomp,ng);
321  MultiFab::Copy(*m_rot_fields[lev][3],*Theta_prim,0,0,1,ng);
322  if (Qv_prim) MultiFab::Copy(*m_rot_fields[lev][4],*Qv_prim,0,0,1,ng);
323  } else {
324  m_rot_fields[lev][3] = nullptr;
325  m_rot_fields[lev][4] = nullptr;
326  }
327 
328  // Default zref to 10 and fill will true values later
329  m_zref[lev] = std::make_unique<MultiFab>(ba2d,dm,1,ng);
330  m_zref[lev]->setVal(zref_default);
331 
332  if (use_terrain_fitted_coords && m_norm_vec && m_interp) {
333  m_x_pos[lev] = std::make_unique<MultiFab>(ba2d,dm,ncomp,ng);
334  m_y_pos[lev] = std::make_unique<MultiFab>(ba2d,dm,ncomp,ng);
335  m_z_pos[lev] = std::make_unique<MultiFab>(ba2d,dm,ncomp,ng);
336  } else if (use_terrain_fitted_coords && m_interp) {
337  m_x_pos[lev] = std::make_unique<MultiFab>(ba2d,dm,ncomp,ng);
338  m_y_pos[lev] = std::make_unique<MultiFab>(ba2d,dm,ncomp,ng);
339  m_z_pos[lev] = std::make_unique<MultiFab>(ba2d,dm,ncomp,ng);
340  } else if (use_terrain_fitted_coords && m_norm_vec) {
341  m_i_indx[lev] = std::make_unique<iMultiFab>(ba2d,dm,incomp,ng);
342  m_j_indx[lev] = std::make_unique<iMultiFab>(ba2d,dm,incomp,ng);
343  m_k_indx[lev] = std::make_unique<iMultiFab>(ba2d,dm,incomp,ng);
344  } else {
345  if (!use_eb) {
346  m_k_indx[lev] = std::make_unique<iMultiFab>(ba2d,dm,incomp,ng);
347  }
348  }
349  }
350 
351  // Setup auxiliary data for spatial configuration & policy
352  //--------------------------------------------------------
353  if (use_terrain_fitted_coords && m_norm_vec && m_interp) { // Terrain w/ norm & w/ interpolation
354  set_norm_positions_T(lev);
355  } else if (use_terrain_fitted_coords && m_interp) { // Terrain w/ interpolation
356  set_z_positions_T(lev);
357  } else if (use_terrain_fitted_coords && m_norm_vec) { // Terrain w/ norm & w/o interpolation
358  set_norm_indices_T(lev);
359  } else if (use_terrain_fitted_coords) { // Terrain
360  set_k_indices_T(lev);
361  } else if (use_eb) { // EB
362  set_z_positions_EB(lev);
363  } else { // No Terrain
364  set_k_indices_N(lev);
365  }
366 
367  // Check once that the direct lateral regional path can read the reference
368  // stencil from every face-owned FAB; distributed staging can replace this
369  // guard in the future.
370  validate_lateral_reference_stencil(lev);
371  // Report the reference height the surface layer will use at this level
372  {
373  const Real zref_min = m_zref[lev]->min(0);
374  const Real zref_max = m_zref[lev]->max(0);
375  Print() << "MOST reference height at level " << lev << ": " << zref_min;
376  if (zref_max > zref_min) { Print() << " to " << zref_max; }
377  Print() << std::endl;
378  }
379 
380  // Setup normalization data for the chosen policy
381  //--------------------------------------------------------
382  switch(m_policy) {
383  case 0: // Plane average
384  set_plane_normalization(lev);
385  break;
386  case 1: // Local region/point
387  set_region_normalization(lev);
388  break;
389  case 2: // EB
390  set_eb_normalization(lev);
391  break;
392  default:
393  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(false, "Unknown policy for MOSTAverage!");
394  }
395 
396  // Set up the exponential time filtering
397  //--------------------------------------------------------
398  if (m_t_avg) {
399  // Exponential filter function
400  m_fact_old = static_cast<amrex::Real>(std::exp(-1.0 / m_time_window));
401 
402  // Enforce discrete normalization: (mfn*val_new + mfo*val_old)
403  m_fact_new = one - m_fact_old;
404 
405  // None of the averages are initialized
406  m_t_init.resize(m_maxlev,0);
407 
408  // We have just (re)built the average containers at this level, so whatever
409  // filter history this level held is gone: m_averages holds bogus_large_value
410  // and set_plane_normalization has zeroed m_plane_average. Note that resize
411  // above does not touch existing entries, so this must be set explicitly or a
412  // regrid would blend the bogus values into the filtered average. On restart
413  // ReadCheckpointFileSurfaceLayer restores the state and sets this back to 1.
414  m_t_init[lev] = 0;
415  }
416 
417  // Correction to the mean surface velocity at this level
418  m_Vsg[lev] = zero;
419  if (include_subgrid_vel) {
420  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(dir == 2, "Subgrid velocity correction only for Z faces!");
421  Print() << "Subgrid velocity scale correction at level : " << lev << ' ';
422  const auto dxArr = m_geom[lev].CellSizeArray();
423  Real dx = std::sqrt(dxArr[0]*dxArr[1]);
424  if (dx > Real(5000.)) {
425  m_Vsg[lev] = Real(0.32) * std::pow(dx/Real(5000.)-1, Real(0.33));
426  }
427  Print() << m_Vsg[lev] << std::endl;
428  }
429 }
bogus_large_value
constexpr amrex::Real bogus_large_value
Definition: ERF_Constants.H:19
use_eb
const bool use_eb
Definition: ERF_InitCustomPert_KE.H:8
MOSTAverage::m_time_window
double m_time_window
Definition: ERF_MOSTAverage.H:606
MOSTAverage::set_z_positions_EB
void set_z_positions_EB(const int &lev)
Definition: ERF_MOSTAverage.cpp:979
MOSTAverage::validate_lateral_reference_stencil
void validate_lateral_reference_stencil(const int &lev)
Definition: ERF_MOSTAverage.cpp:940
MOSTAverage::set_region_normalization
void set_region_normalization(const int &)
Definition: ERF_MOSTAverage.H:131
MOSTAverage::set_z_positions_T
void set_z_positions_T(const int &lev)
Definition: ERF_MOSTAverage.cpp:1275
MOSTAverage::set_norm_positions_T
void set_norm_positions_T(const int &lev)
Definition: ERF_MOSTAverage.cpp:1354
MOSTAverage::set_k_indices_T
void set_k_indices_T(const int &lev)
Definition: ERF_MOSTAverage.cpp:1007
MOSTAverage::set_plane_normalization
void set_plane_normalization(const int &lev)
Definition: ERF_MOSTAverage.cpp:581
MOSTAverage::set_norm_indices_T
void set_norm_indices_T(const int &lev)
Definition: ERF_MOSTAverage.cpp:1162
MOSTAverage::m_maxlev
int m_maxlev
Definition: ERF_MOSTAverage.H:567
MOSTAverage::set_eb_normalization
void set_eb_normalization(const int &lev)
Definition: ERF_MOSTAverage.cpp:639
MOSTAverage::include_subgrid_vel
bool include_subgrid_vel
Definition: ERF_MOSTAverage.H:611
MOSTAverage::set_k_indices_N
void set_k_indices_N(const int &lev)
Definition: ERF_MOSTAverage.cpp:733
MOSTAverage::zref_default
const amrex::Real zref_default
Definition: ERF_MOSTAverage.H:616
Vars::xvel
@ xvel
Definition: ERF_IndexDefines.H:255
Vars::zvel
@ zvel
Definition: ERF_IndexDefines.H:257
Vars::yvel
@ yvel
Definition: ERF_IndexDefines.H:256

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

Deleted copy-assignment operator.

Parameters
[in]othersource object

◆ operator=() [2/2]

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

Deleted move-assignment operator.

Parameters
[in]othersource object

◆ set_eb_normalization()

void MOSTAverage::set_eb_normalization ( const int &  lev)

Compute total embedded-boundary surface area.

Parameters
[in]levlevel index

Function to compute normalization for average over EB.

Parameters
[in]levCurrent level
640 {
641  AMREX_ALWAYS_ASSERT(m_eb_vec[lev] != nullptr);
642 
643  // Get EB data - need both cell-centered and face-centered flags
644  const auto& cc_flags = m_eb_vec[lev]->get_const_factory()->getMultiEBCellFlagFab();
645 
646  // Get area fractions for different centerings
647  auto cc_afrac = m_eb_vec[lev]->get_const_factory()->getAreaFrac(); // Cell-centered area fractions (Array of 3 MultiCutFab*)
648 
649  // Initialize storage
650  m_total_bndry_area.resize(m_maxlev);
651  m_total_bndry_area[lev].resize(m_navg, zero);
652  m_plane_average.resize(m_maxlev);
653  m_plane_average[lev].resize(m_navg, zero);
654 
655  // Compute total area for each field type based on its centering
656  // iavg: 0=U(xface), 1=V(yface), 2=W(zface), 3=T(cc), 4=Qv(cc),
657  // 5=Tv(cc), 6=Umag_XY(cc), 7=Umag_XZ(cc), 8=Umag_YZ(cc)
658 
659  // Get geometry for cell sizes
660  auto const& dx_arr = m_geom[lev].CellSizeArray();
661  Real dx = dx_arr[0];
662  Real dy = dx_arr[1];
663  Real dz = dx_arr[2];
664 
665  // GPU array to accumulate areas
666  Gpu::DeviceVector<Real> area_vec(m_navg, zero);
667  Real* area_device = area_vec.data();
668 
669  // All fields are now averaged on cell-centered grid
670  // Compute total area once on cell-centered grid and use for all iavg
671  Real total_area = zero;
672  //
673  // NOTE: deliberately not threaded. The body is a pure reduction into the shared
674  // area_device entry, and on the host Gpu::deviceReduceSum is an unordered
675  // "#pragma omp atomic" add. Threading it would make the cut-cell area -- which
676  // normalizes every EB surface-layer average -- depend on thread scheduling in the
677  // last bits. See the same note in PlaneAverage::compute_averages.
678  //
679  for (MFIter mfi(cc_flags, TileNoZ()); mfi.isValid(); ++mfi) {
680  const auto& flag = cc_flags[mfi];
681 
682  // Skip boxes that are not singlevalued (MultiCutFab only has data for singlevalued boxes)
683  if (flag.getType() != FabType::singlevalued) continue;
684 
685  Box bx = mfi.tilebox();
686  auto const flag_arr = flag.const_array();
687  auto const afrac_x = cc_afrac[0]->const_array(mfi);
688  auto const afrac_y = cc_afrac[1]->const_array(mfi);
689  auto const afrac_z = cc_afrac[2]->const_array(mfi);
690 
691  // Sum area only for cut cells using atomic reduction
692  ParallelFor(Gpu::KernelInfo().setReduction(true), bx, [=]
693  AMREX_GPU_DEVICE(int i, int j, int k, Gpu::Handler const& handler) noexcept
694  {
695  if (flag_arr(i,j,k).isSingleValued()) {
696  // Compute area from face-centered area fractions
697  Real axm = afrac_x(i ,j ,k );
698  Real axp = afrac_x(i+1,j ,k );
699  Real aym = afrac_y(i ,j ,k );
700  Real ayp = afrac_y(i ,j+1,k );
701  Real azm = afrac_z(i ,j ,k );
702  Real azp = afrac_z(i ,j ,k+1);
703 
704  Real adx = (axm - axp) * dy * dz;
705  Real ady = (aym - ayp) * dx * dz;
706  Real adz = (azm - azp) * dx * dy;
707 
708  Real area = std::sqrt(adx*adx + ady*ady + adz*adz);
709 
710  Gpu::deviceReduceSum(&area_device[0], area, handler);
711  }
712  });
713  }
714 
715  // Copy to host and sum across MPI ranks
716  Gpu::copy(Gpu::deviceToHost, area_vec.begin(), area_vec.begin() + 1, &total_area);
717  ParallelDescriptor::ReduceRealSum(&total_area, 1);
718 
719  // Set the same total area for all iavg
720  for (int iavg = 0; iavg < m_navg; ++iavg) {
721  m_total_bndry_area[lev][iavg] = total_area;
722  }
723 
724  Print() << "EB surface area on cell-centerd grid at level " << lev << ": " << total_area << std::endl;
725 }

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)

Populate the 2D k-index iMultiFab without terrain.

Parameters
[in]levlevel index

Function to set K indices without terrain.

Parameters
[in]levCurrent level
734 {
735  ParmParse pp(m_pp_prefix);
736  // See zref_sentinel in ERF_MOSTAverage.H: probe with the sentinel and test the
737  // value, since queryAdd's return value stops meaning "user specified" after the
738  // first parse of the key and this routine runs once per level.
739  Real zref_tmp = zref_sentinel;
740  pp.queryAdd("most.zref",zref_tmp);
741  bool read_z = (zref_tmp > Real(0));
742  if (!read_z) { zref_tmp = zref_default; }
743  auto read_k = pp.queryarr("most.k_arr_in",m_k_in);
744 
745  const int dir = m_face.coordDir();
746  const bool is_lo_face = m_face.isLow();
747  const bool zlo = (dir == 2 && is_lo_face);
748 
749  const bool stretched = (dir == 2 && m_mesh_type != MeshType::ConstantDz);
750  if (stretched) {
751  for (int ilev : {0, lev}) {
752  const int nz = m_geom[ilev].Domain().length(2);
753  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
754  ilev < static_cast<int>(m_zlevels_stag.size()) &&
755  static_cast<int>(m_zlevels_stag[ilev].size()) == nz + 1,
756  "MOSTAverage: z levels are needed for a stretched mesh without terrain!");
757  }
758  }
759 
760  // Height above the wall of the default reference on zlo of a constant-dz mesh, kept
761  // apart from the absolute z the input convention there uses, so that subtracting
762  // prob_lo again cannot round it below the first cell center (negative when unset)
763  Real zlo_default_height = -one;
764 
765  // Default behavior is to use the first cell center
766  if (!read_z && !read_k) {
767  if (stretched) {
768  const auto& zlevels = m_zlevels_stag[0];
769  zref_tmp = zlo ? cell_center_height(zlevels, 0)
770  : column_height(zlevels) - cell_center_height(zlevels, zlevels.size()-2);
771  } else {
772  const Real dz = m_geom[0].CellSize(dir);
773  zref_tmp = zlo ? m_geom[0].ProbLo(dir) + myhalf * dz : myhalf * dz;
774  if (zlo) { zlo_default_height = myhalf * dz; }
775  }
776  if (zlo && !stretched) {
777  // zlo input convention: an absolute z
778  Print() << "Reference height for MOST set to " << zref_tmp - m_geom[0].ProbLo(dir)
779  << " (z = " << zref_tmp << ")" << std::endl;
780  } else {
781  Print() << "Reference height for MOST set to " << zref_tmp << std::endl;
782  }
783  read_z = true;
784  }
785 
786  // Specify z_ref & compute k_indx (z_ref takes precedence)
787  if (read_z && stretched) {
788  const auto& zlevels = m_zlevels_stag[lev];
789  const int nz = static_cast<int>(zlevels.size()) - 1;
790  const Real ztop = column_height(zlevels);
791  const Real target = zlo ? zref_tmp : ztop - zref_tmp;
792  const Real first_center = cell_center_height(zlevels, 0);
793  const Real last_center = cell_center_height(zlevels, nz - 1);
794  const Real target_min = zlo ? first_center : ztop - last_center;
795  const Real target_max = zlo ? last_center : ztop - first_center;
796  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(zref_tmp >= target_min && zref_tmp <= target_max,
797  "Query point must remain inside the stretched mesh!");
798  const int lk = k_index_below(zlevels, target);
799  AMREX_ALWAYS_ASSERT_WITH_MESSAGE((zlo ? lk : nz - 1 - lk) >= m_radius,
800  "K index must be larger than averaging radius!");
801  m_k_indx[lev]->setVal(lk);
802  m_zref[lev]->setVal(zlo ? cell_center_height(zlevels, lk)
803  : ztop - cell_center_height(zlevels, lk));
804  } else if (read_z) {
805  const Real m_zlo = m_geom[lev].ProbLo(dir);
806  const Real m_zhi = m_geom[lev].ProbHi(dir);
807  const Real m_dz = m_geom[lev].CellSize(dir);
808 
809  const int dom_lo = m_geom[lev].Domain().smallEnd(dir);
810  const int dom_hi = m_geom[lev].Domain().bigEnd(dir);
811 
812  if (zlo) {
813  // Existing behavior: zlo input is an absolute coordinate. The stored
814  // reference height is the height above the wall, as on every other face:
815  // the surface layer uses it as one, e.g. in ln(zref / z0).
816  // Height above the wall. A query at a cell center (prob_lo + dz/2 typed as a z,
817  // or the default dz0/2 on a level refined by an odd ratio) must not round below
818  // it, so the cell count gets an allowance of a few round-off units of the
819  // numbers it came from
820  const Real height = (zlo_default_height >= zero) ? zlo_default_height
821  : zref_tmp - m_zlo;
822  const Real slack = Real(8) * std::numeric_limits<Real>::epsilon()
823  * ((std::abs(zref_tmp) + std::abs(m_zlo)) / m_dz + one);
824  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(height / m_dz + slack >= myhalf,
825  "erf.most.zref on the low z face of a constant-dz mesh is an absolute z: it must be "
826  "at least geometry.prob_lo[2] + dz/2");
827  AMREX_ALWAYS_ASSERT_WITH_MESSAGE((m_zhi - m_zlo - height) / m_dz + slack >= myhalf,
828  "erf.most.zref on the low z face of a constant-dz mesh is an absolute z: it must be "
829  "at most geometry.prob_hi[2] - dz/2");
830 
831  const int wall_offset =
832  static_cast<int>(floor(height / m_dz - myhalf + slack));
833 
834  AMREX_ALWAYS_ASSERT(wall_offset >= m_radius);
835 
836  m_k_indx[lev]->setVal(dom_lo + wall_offset);
837  m_zref[lev]->setVal((static_cast<Real>(wall_offset) + myhalf) * m_dz);
838  } else {
839  const Real domain_length = m_zhi - m_zlo;
840 
841  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(zref_tmp >= myhalf * m_dz,
842  "Query distance must reach the first cell center!");
843  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(zref_tmp <= domain_length - myhalf * m_dz,
844  "Query distance must remain inside the domain!");
845 
846  const int wall_offset = static_cast<int>(floor(zref_tmp / m_dz - myhalf));
847  const int lk = is_lo_face ? dom_lo + wall_offset : dom_hi - wall_offset;
848 
849  AMREX_ALWAYS_ASSERT(wall_offset >= m_radius);
850 
851  const Real zref_abs = is_lo_face
852  ? m_zlo + (static_cast<Real>(wall_offset) + myhalf) * m_dz
853  : m_zhi - (static_cast<Real>(wall_offset) + myhalf) * m_dz;
854  const Real zref_height = is_lo_face
855  ? zref_abs - m_zlo
856  : m_zhi - zref_abs;
857 
858  m_k_indx[lev]->setVal(lk);
859  m_zref[lev]->setVal(zref_height);
860  }
861  // Specified k_indx & compute z_ref
862  } else if (read_k) {
863  const int dom_lo = m_geom[lev].Domain().smallEnd(dir);
864  const int dom_hi = m_geom[lev].Domain().bigEnd(dir);
865  const int ncell = m_geom[lev].Domain().length(dir);
866 
867  // Preserve the established z-low convention: k_arr_in is an
868  // absolute cell index on zlo. All other faces interpret k_arr_in as
869  // a distance in cells from the selected wall.
870  const int wall_offset = zlo
871  ? m_k_in[lev] - dom_lo
872  : m_k_in[lev];
873  if (zlo) {
874  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
875  m_k_in[lev] >= dom_lo && m_k_in[lev] <= dom_hi,
876  "MOST zlo reference index must lie inside the domain!");
877  } else {
878  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
879  wall_offset >= 0 && wall_offset < ncell,
880  "MOST wall-relative reference offset must lie inside the domain!");
881  }
882  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(wall_offset >= m_radius,
883  "K index must be larger than averaging radius!");
884  const int ref_index = zlo
885  ? m_k_in[lev]
886  : (is_lo_face ? dom_lo + wall_offset : dom_hi - wall_offset);
887  m_k_indx[lev]->setVal(ref_index);
888 
889  // TODO: check that z_ref is constant across levels
890  if (dir == 2) {
891  const int dom_lo0 = m_geom[0].Domain().smallEnd(2);
892  const int dom_hi0 = m_geom[0].Domain().bigEnd(2);
893  const int wall_offset0 = zlo ? m_k_in[0] - dom_lo0 : m_k_in[0];
894  const int ref_index0 = zlo ? m_k_in[0] : dom_hi0 - wall_offset0;
895  const int k0 = ref_index0 - dom_lo0;
896  AMREX_ALWAYS_ASSERT(wall_offset0 >= m_radius);
897 
898  if (stretched) {
899  const auto& zlevels = m_zlevels_stag[0];
900  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
901  k0 >= 0 && k0 < m_geom[0].Domain().length(2),
902  "MOST reference index must lie inside the stretched domain!");
903  const Real zcell = cell_center_height(zlevels, k0);
904  m_zref[lev]->setVal(zlo ? zcell : column_height(zlevels) - zcell);
905  } else {
906  const Real dz0 = m_geom[0].CellSize(2);
907  const Real zlo0 = m_geom[0].ProbLo(2);
908  const Real zhi0 = m_geom[0].ProbHi(2);
909  const Real zref_abs = is_lo_face
910  ? zlo0 + (static_cast<Real>(wall_offset0) + myhalf) * dz0
911  : zhi0 - (static_cast<Real>(wall_offset0) + myhalf) * dz0;
912  m_zref[lev]->setVal(is_lo_face ? zref_abs - zlo0 : zhi0 - zref_abs);
913  }
914  } else {
915  const Real m_dz = m_geom[lev].CellSize(dir);
916  const Real m_zlo = m_geom[lev].ProbLo(dir);
917  const Real m_zhi = m_geom[lev].ProbHi(dir);
918  const Real zref_abs = is_lo_face
919  ? m_zlo + (static_cast<Real>(wall_offset) + myhalf) * m_dz
920  : m_zhi - (static_cast<Real>(wall_offset) + myhalf) * m_dz;
921 
922  m_zref[lev]->setVal(is_lo_face ? zref_abs - m_zlo : m_zhi - zref_abs);
923  }
924  }
925 }
ztop
const Real ztop
Definition: ERF_InitCustomPertVels_ParticleTests.H:4
length
const Real length
Definition: ERF_InitCustomPert_AnelasticWallDiffusion.H:16
pp
ParmParse pp("prob")
height
Real height
Definition: ERF_InitCustomPert_SquallLine.H:37
MOSTAverage::column_height
static amrex::Real column_height(const amrex::Vector< amrex::Real > &zlevels)
Definition: ERF_MOSTAverage.H:396
MOSTAverage::k_index_below
static int k_index_below(const amrex::Vector< amrex::Real > &zlevels, const amrex::Real &zref)
Definition: ERF_MOSTAverage.H:366
MOSTAverage::zref_sentinel
const amrex::Real zref_sentinel
Definition: ERF_MOSTAverage.H:626
MOSTAverage::m_pp_prefix
std::string m_pp_prefix
Definition: ERF_MOSTAverage.H:558
MOSTAverage::m_mesh_type
MeshType m_mesh_type
Definition: ERF_MOSTAverage.H:559
MOSTAverage::m_zlevels_stag
amrex::Vector< amrex::Vector< amrex::Real > > m_zlevels_stag
Definition: ERF_MOSTAverage.H:561
MOSTAverage::m_k_in
amrex::Vector< int > m_k_in
Definition: ERF_MOSTAverage.H:590
module_model_constants::epsilon
real(c_double), parameter epsilon
Definition: ERF_module_model_constants.F90:12

Referenced by make_MOSTAverage_at_level(), and set_k_indices_T().

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)

Populate the 2D k-index iMultiFab with terrain.

Parameters
[in]levlevel index

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

Parameters
[in]levCurrent level
1008 {
1009  // Terrain-fitted coordinates only deform the vertical coordinate. A
1010  // lateral wall therefore uses the uniform wall-normal indexing logic.
1011  if (m_face.coordDir() != 2) {
1012  set_k_indices_N(lev);
1013  return;
1014  }
1015 
1016  ParmParse pp(m_pp_prefix);
1017  // See zref_sentinel in ERF_MOSTAverage.H: probe with the sentinel and test the
1018  // value, since queryAdd's return value stops meaning "user specified" after the
1019  // first parse of the key and this routine runs once per level.
1020  Real zref_tmp = zref_sentinel;
1021  pp.queryAdd("most.zref",zref_tmp);
1022  const bool user_zref = (zref_tmp > Real(0));
1023  bool read_z = user_zref;
1024  auto read_k = pp.queryarr("most.k_arr_in",m_k_in);
1025 
1026  // Allow default zref
1027  if (!read_z) {
1028  zref_tmp = zref_default;
1029  Print() << "most.zref not specified, query distance default is " << zref_tmp << std::endl;
1030  read_z = true;
1031  }
1032 
1033  // No default behavior with terrain (we can't tell the difference between
1034  // vertical grid stretching and true terrain)
1035  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(read_z != read_k,
1036  "Need to specify zref or k_arr_in for MOST");
1037 
1038  // Capture for device
1039  Real d_zref = zref_tmp;
1040  Real d_radius = static_cast<Real>(m_radius);
1041  const bool is_lo_face = m_face.isLow();
1042  const int ng = m_k_indx[lev]->nGrow() - 1;
1043 
1044  // The k indices are needed everywhere an average is computed, ghost cells
1045  // included; the box is made nodal so that we also cover the U & V averages,
1046  // which are face centered
1047  const IntVect ng_indx = max(get_ng_fill(lev), IntVect(ng,ng,0));
1048 
1049  // Specify z_ref & compute k_indx (z_ref takes precedence)
1050  if (read_z) {
1051  const Box domain = m_geom[lev].Domain();
1052  const int zlo = domain.smallEnd(2);
1053  const int zhi = domain.bigEnd(2);
1054  const int top_node = zhi + 1;
1055 
1056  // Iterate the full 3-D source layout so z-split grids can be tested
1057  // for ownership of the selected wall before collapsing the kernel box.
1058  for (MFIter mfi(*m_fields[lev][3], TileNoZ()); mfi.isValid(); ++mfi) {
1059  Box npbx = mfi.tilebox(IntVect(1,1,0),ng_indx);
1060  const Box vbx = mfi.validbox();
1061 
1062  if (is_lo_face) {
1063  if (vbx.smallEnd(2) != zlo) { continue; }
1064  npbx.makeSlab(2,zlo);
1065  } else {
1066  if (vbx.bigEnd(2) != zhi) { continue; }
1067  npbx.makeSlab(2,zhi);
1068  }
1069 
1070  const Box& z_phys_fab = m_z_phys_nd[lev]->fabbox(mfi.index());
1071  // The z_phys FAB may cover only one slab of a vertically
1072  // decomposed domain. The zlo search starts at the bottom and
1073  // normally finds the reference cell before leaving that FAB, so
1074  // it must not require the domain-wide top node. The zhi search
1075  // does require the top node because that is its first interval.
1076  if (!is_lo_face) {
1077  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
1078  z_phys_fab.bigEnd(2) >= top_node,
1079  "Terrain FAB does not contain the top nodal plane.");
1080  }
1081  const int search_last_cell = is_lo_face
1082  ? min(zhi, z_phys_fab.bigEnd(2) - 1)
1083  : max(zlo, z_phys_fab.smallEnd(2));
1084  const int first_search_cell = is_lo_face ? zlo : zhi;
1085  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
1086  (is_lo_face && search_last_cell >= first_search_cell) ||
1087  (!is_lo_face && search_last_cell <= first_search_cell),
1088  "Terrain FAB does not contain a complete search interval.");
1089  const auto z_phys_arr = m_z_phys_nd[lev]->const_array(mfi);
1090  auto k_arr = m_k_indx[lev]->array(mfi);
1091  auto zref_arr = m_zref[lev]->array(mfi);
1092  auto z_at = [=] AMREX_GPU_DEVICE (int i, int j, int node_k) noexcept -> Real
1093  {
1094  return fourth * ( z_phys_arr(i ,j ,node_k) + z_phys_arr(i+1,j ,node_k)
1095  + z_phys_arr(i ,j+1,node_k) + z_phys_arr(i+1,j+1,node_k) );
1096  };
1097 
1098  ParallelFor(npbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
1099  {
1100  // Terrain heights are nodal. Use the horizontal average of
1101  // the four nodes at a given vertical node for the face and
1102  // for every cell interval searched below.
1103  const int face_node = is_lo_face ? zlo : top_node;
1104  const Real z_face = z_at(i, j, face_node);
1105  const Real z_target = z_face + (is_lo_face ? d_zref : -d_zref);
1106  const int first_cell = is_lo_face ? zlo : zhi;
1107  const int cell_step = is_lo_face ? 1 : -1;
1108  k_arr(i,j,k) = first_cell;
1109 
1110  bool found = false;
1111  for (int cell = first_cell; ; cell += cell_step) {
1112  const Real z_lo = z_at(i, j, cell);
1113  const Real z_hi = z_at(i, j, cell + 1);
1114  // Preserve the legacy wall-adjacent choice for an explicitly
1115  // requested zref. The half-open interval is needed for the
1116  // default query on stretched meshes, where the default can
1117  // land exactly on a cell face.
1118  const bool in_cell = user_zref
1119  ? (z_target >= z_lo && z_target <= z_hi)
1120  : (is_lo_face
1121  ? (z_target >= z_lo && z_target < z_hi)
1122  : (z_target > z_lo && z_target <= z_hi));
1123  if (in_cell) {
1124  const int wall_offset = is_lo_face ? cell - zlo : zhi - cell;
1125  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(wall_offset >= d_radius,
1126  "K index must be larger than averaging radius!");
1127  k_arr(i,j,k) = cell;
1128  const Real z_cell = myhalf * (z_hi + z_lo);
1129  zref_arr(i,j,k) = is_lo_face ? z_cell - z_face : z_face - z_cell;
1130  found = true;
1131  break;
1132  }
1133  // z_at(cell) reads both cell nodes, so stop at the last
1134  // interval available in this FAB. If the reference
1135  // height is not available locally, the assertion below
1136  // reports the problem instead of reading another slab.
1137  if (cell == search_last_cell) { break; }
1138  }
1139 
1140  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(found,
1141  "zref not found with terrain!");
1142  });
1143  }
1144 
1145  const Periodicity period = tangential_periodicity(m_geom[lev], 2);
1146  m_k_indx[lev]->FillBoundary(m_k_indx[lev]->nGrowVect(), period);
1147  m_zref[lev]->FillBoundary(m_zref[lev]->nGrowVect(), period);
1148 
1149  // Specified k_indx & compute z_ref
1150  } else if (read_k) {
1151  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(false, "Specified k-indx with terrain not implemented!");
1152  }
1153 }
fourth
constexpr amrex::Real fourth
Definition: ERF_NumericalConstants.H:42

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)

Populate all 2D normal-index iMultiFabs with terrain.

Parameters
[in]levlevel index

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

Parameters
[in]levCurrent level
1163 {
1164  // Peel back the level
1165  auto& fields = m_fields[lev];
1166 
1167  // MFIter over CC data
1168  int imf_cc = 3;
1169 
1170  ParmParse pp(m_pp_prefix);
1171  // See zref_sentinel in ERF_MOSTAverage.H: probe with the sentinel and test the
1172  // value, since queryAdd's return value stops meaning "user specified" after the
1173  // first parse of the key and this routine runs once per level.
1174  Real zref_tmp = zref_sentinel;
1175  pp.queryAdd("most.zref",zref_tmp);
1176  bool read_zref = (zref_tmp > Real(0));
1177  if (!read_zref) {
1178  zref_tmp = zref_default;
1179  Print() << "most.zref not specified, query distance default is " << zref_tmp << std::endl;
1180  }
1181 
1182  AMREX_ALWAYS_ASSERT(m_face.coordDir() == 2);
1183  int klo = m_geom[lev].Domain().smallEnd(2);
1184 
1185  // Capture for device
1186  Real d_zref = zref_tmp;
1187  Real d_radius = static_cast<Real>(m_radius);
1188 
1189  const auto dxInv = m_geom[lev].InvCellSizeArray();
1190 
1191  // The indices are needed everywhere an average is computed, ghost cells
1192  // included; the box is made nodal so that we also cover the U & V averages,
1193  // which are face centered
1194  const IntVect ng_indx = max(get_ng_fill(lev), IntVect(1,1,0));
1195 
1196  for (MFIter mfi(*fields[imf_cc], TileNoZ()); mfi.isValid(); ++mfi) {
1197  Box npbx = mfi.tilebox(IntVect(1,1,0),IntVect(1,1,0));
1198  const Box vbx = mfi.validbox();
1199  npbx = mfi.tilebox(IntVect(1,1,0),ng_indx);
1200 
1201  if (vbx.smallEnd(2) != klo) { continue; }
1202 
1203  int kmax = npbx.bigEnd(2);
1204 
1205  npbx.makeSlab(2,klo);
1206 
1207  const auto z_phys_arr = m_z_phys_nd[lev]->const_array(mfi);
1208  auto i_arr = m_i_indx[lev]->array(mfi);
1209  auto j_arr = m_j_indx[lev]->array(mfi);
1210  auto k_arr = m_k_indx[lev]->array(mfi);
1211  auto zref_arr = m_zref[lev]->array(mfi);
1212 
1213  // The displaced column (i_new, j_new) is read here from z_phys at i_new..i_new+1 and
1214  // j_new..j_new+1, and later from the averaged fields at (i_new, j_new); keep it
1215  // inside both FABs, so a displacement past the ghost cells uses the nearest column
1216  const auto zp_lo = lbound(z_phys_arr), zp_hi = ubound(z_phys_arr);
1217  const auto cc_lo = lbound(fields[imf_cc]->const_array(mfi));
1218  const auto cc_hi = ubound(fields[imf_cc]->const_array(mfi));
1219  const int inew_lo = amrex::max(zp_lo.x, cc_lo.x), inew_hi = amrex::min(zp_hi.x - 1, cc_hi.x);
1220  const int jnew_lo = amrex::max(zp_lo.y, cc_lo.y), jnew_hi = amrex::min(zp_hi.y - 1, cc_hi.y);
1221 
1222  ParallelFor(npbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
1223  {
1224  // Elements of normal vector
1225  Real met_h_xi = Compute_h_xi_AtCellCenter (i,j,k,dxInv,z_phys_arr);
1226  Real met_h_eta = Compute_h_eta_AtCellCenter(i,j,k,dxInv,z_phys_arr);
1227  Real mag = std::sqrt(met_h_xi*met_h_xi + met_h_eta*met_h_eta + one);
1228 
1229  // Unit-normal vector scaled by z_ref
1230  Real delta_x = -met_h_xi/mag * d_zref;
1231  Real delta_y = -met_h_eta/mag * d_zref;
1232  Real delta_z = one/mag * d_zref;
1233 
1234  // Compute i & j as displacements (no grid stretching)
1235  int delta_i = static_cast<int>(std::round(delta_x*dxInv[0]));
1236  int delta_j = static_cast<int>(std::round(delta_y*dxInv[1]));
1237  int i_new = amrex::Clamp(i + delta_i, inew_lo, inew_hi);
1238  int j_new = amrex::Clamp(j + delta_j, jnew_lo, jnew_hi);
1239  i_arr(i,j,0) = i_new;
1240  j_arr(i,j,0) = j_new;
1241 
1242  // Search for k (grid is stretched in z)
1243  Real z_bot_face = fourth * ( z_phys_arr(i ,j ,k) + z_phys_arr(i+1,j ,k)
1244  + z_phys_arr(i ,j+1,k) + z_phys_arr(i+1,j+1,k) );
1245  Real z_target = z_bot_face + delta_z;
1246  k_arr(i,j,0) = klo;
1247  zref_arr(i,j,0) = myhalf * z_bot_face +
1248  Real(0.125) * ( z_phys_arr(i ,j ,k+1) + z_phys_arr(i+1,j ,k+1)
1249  + z_phys_arr(i ,j+1,k+1) + z_phys_arr(i+1,j+1,k+1) );
1250  for (int lk(klo); lk<=kmax; ++lk) {
1251  Real z_lo = fourth * ( z_phys_arr(i_new,j_new ,lk ) + z_phys_arr(i_new+1,j_new ,lk )
1252  + z_phys_arr(i_new,j_new+1,lk ) + z_phys_arr(i_new+1,j_new+1,lk ) );
1253  Real z_hi = fourth * ( z_phys_arr(i_new,j_new ,lk+1) + z_phys_arr(i_new+1,j_new ,lk+1)
1254  + z_phys_arr(i_new,j_new+1,lk+1) + z_phys_arr(i_new+1,j_new+1,lk+1) );
1255  if (in_cell_z(z_target, z_lo, z_hi)) {
1256  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(lk >= d_radius,
1257  "K index must be larger than averaging radius!");
1258  amrex::ignore_unused(d_radius);
1259  k_arr(i,j,0) = lk;
1260  zref_arr(i,j,0) = myhalf * (z_hi + z_lo) - z_bot_face;
1261  break;
1262  }
1263  }
1264  });
1265  }
1266 }
klo
const int klo
Definition: ERF_InitCustomPert_KE.H:9
Compute_h_xi_AtCellCenter
AMREX_GPU_HOST_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:226
Compute_h_eta_AtCellCenter
AMREX_GPU_HOST_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:251
MOSTAverage::in_cell_z
AMREX_GPU_HOST_DEVICE static AMREX_INLINE bool in_cell_z(const amrex::Real &z, const amrex::Real &z_lo, const amrex::Real &z_hi)
Definition: ERF_MOSTAverage.H:350

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)

Populate terrain-aware normal positions.

Parameters
[in]levlevel index

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

Parameters
[in]levCurrent level
1355 {
1356  AMREX_ALWAYS_ASSERT(m_face.coordDir() == 2);
1357 
1358  // Peel back the level
1359  auto& fields = m_fields[lev];
1360 
1361  // MFIter over CC data
1362  int imf_cc = 3;
1363 
1364  ParmParse pp(m_pp_prefix);
1365  // See zref_sentinel in ERF_MOSTAverage.H: probe with the sentinel and test the
1366  // value, since queryAdd's return value stops meaning "user specified" after the
1367  // first parse of the key and this routine runs once per level.
1368  Real zref_tmp = zref_sentinel;
1369  pp.queryAdd("most.zref",zref_tmp);
1370  bool read_zref = (zref_tmp > Real(0));
1371  if (!read_zref) {
1372  zref_tmp = zref_default;
1373  Print() << "most.zref not specified, query distance default is " << zref_tmp << std::endl;
1374  }
1375  m_zref[lev]->setVal(zref_tmp);
1376  int klo = m_geom[lev].Domain().smallEnd(2);
1377 
1378  // Capture for device
1379  Real d_zref = zref_tmp;
1380  const auto plo = m_geom[lev].ProbLoArray();
1381 
1382  RealVect base;
1383  const auto dx = m_geom[lev].CellSizeArray();
1384  const auto dxInv = m_geom[lev].InvCellSizeArray();
1385  IntVect ng = m_x_pos[lev]->nGrowVect(); ng[2]=0;
1386  const int position_ng = (m_radius > 1) ? m_radius : 1;
1387 
1388  // The positions are read over the averaging stencil at every cell where an
1389  // average is computed, ghost cells included; the box is made nodal so that
1390  // we also cover the U & V averages, which are face centered
1391  const IntVect ng_pos = max(get_ng_fill(lev) + IntVect(m_radius,m_radius,0),
1392  IntVect(position_ng,position_ng,0));
1393 
1394  for (MFIter mfi(*fields[imf_cc], TileNoZ()); mfi.isValid(); ++mfi) {
1395  Box npbx = mfi.tilebox(IntVect(1,1,0),ng_pos);
1396  Box gtbx = mfi.growntilebox(ng);
1397  RealBox grb{gtbx,dx.data(),base.dataPtr()};
1398 
1399  if (mfi.validbox().smallEnd(2) != klo) { continue; }
1400 
1401  npbx.makeSlab(2,klo);
1402 
1403  const auto z_phys_arr = m_z_phys_nd[lev]->const_array(mfi);
1404  auto x_pos_arr = m_x_pos[lev]->array(mfi);
1405  auto y_pos_arr = m_y_pos[lev]->array(mfi);
1406  auto z_pos_arr = m_z_pos[lev]->array(mfi);
1407  auto zref_arr = m_zref[lev]->array(mfi);
1408  ParallelFor(npbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
1409  {
1410  // Elements of normal vector
1411  Real met_h_xi = Compute_h_xi_AtCellCenter (i,j,k,dxInv,z_phys_arr);
1412  Real met_h_eta = Compute_h_eta_AtCellCenter(i,j,k,dxInv,z_phys_arr);
1413  Real imag = one / std::sqrt(met_h_xi*met_h_xi + met_h_eta*met_h_eta + one);
1414 
1415  // Unit-normal vector scaled by z_ref
1416  Real delta_x = -met_h_xi * imag * d_zref;
1417  Real delta_y = -met_h_eta * imag * d_zref;
1418  Real delta_z = imag * d_zref;
1419 
1420  // Position of the current node (indx:0,0,1)
1421  Real x0 = plo[0] + ((Real) i + myhalf) * dx[0];
1422  Real y0 = plo[1] + ((Real) j + myhalf) * dx[1];
1423 
1424  // Final position at end of vector
1425  x_pos_arr(i,j,0) = x0 + delta_x;
1426  y_pos_arr(i,j,0) = y0 + delta_y;
1427  Real z_bot_face = fourth * ( z_phys_arr(i ,j ,k) + z_phys_arr(i+1,j ,k)
1428  + z_phys_arr(i ,j+1,k) + z_phys_arr(i+1,j+1,k) );
1429  z_pos_arr(i,j,0) = z_bot_face + delta_z;
1430 
1431  // NOTE: Normal vector end point can be below the surface for concave regions.
1432  // Here we protect against that by augmenting the normal if needed.
1433  int i_new = (int) ((x_pos_arr(i,j,0) - plo[0]) / dx[0] - myhalf);
1434  int j_new = (int) ((y_pos_arr(i,j,0) - plo[1]) / dx[1] - myhalf);
1435  Real z_new_bot_face = fourth * ( z_phys_arr(i_new,j_new ,k) + z_phys_arr(i_new+1,j_new ,k)
1436  + z_phys_arr(i_new,j_new+1,k) + z_phys_arr(i_new+1,j_new+1,k) );
1437  if (z_pos_arr(i,j,0) < z_new_bot_face) {
1438  z_pos_arr(i,j,0) = z_new_bot_face + delta_z;
1439  }
1440 
1441  zref_arr(i,j,0) = delta_z;
1442 
1443  // Destination position must be contained on the current process!
1444  Real pos[] = {x_pos_arr(i,j,0)-plo[0],y_pos_arr(i,j,0)-plo[1],myhalf*dx[2]};
1445  amrex::ignore_unused(pos);
1446  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(grb.contains(&pos[0]),
1447  "Query point outside of proc domain!");
1448  });
1449  }
1450 }

Referenced by make_MOSTAverage_at_level().

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

◆ set_plane_average()

bool MOSTAverage::set_plane_average ( const int &  lev,
const amrex::Vector< amrex::Real > &  pavg 
)
inline

Restore the filtered plane averages from a checkpoint. Returns false if the checkpoint does not hold what this run expects, e.g. because the averaging policy changed, in which case nothing is restored.

Parameters
[in]levlevel index
[in]pavgfiltered plane averages, one per average component
316  {
317  const int n_have = (lev < static_cast<int>(m_plane_average.size())) ?
318  static_cast<int>(m_plane_average[lev].size()) : 0;
319  if (n_have != static_cast<int>(pavg.size())) { return false; }
320  if (n_have > 0) { m_plane_average[lev] = pavg; }
321  return true;
322  }

Referenced by SurfaceLayer::set_mac_plane_avg().

Here is the caller graph for this function:

◆ set_plane_normalization()

void MOSTAverage::set_plane_normalization ( const int &  lev)

Compute number of cells per averaging plane.

Parameters
[in]levlevel index

Function to compute normalization for plane average.

Parameters
[in]levCurrent level
582 {
583  // Cells per plane and temp avg storage
584  m_ncell_plane.resize(m_maxlev);
585  m_plane_average.resize(m_maxlev);
586 
587  // True domain not used for normalization
588  Box domain = m_geom[lev].Domain();
589 
590  // NOTE: Level 0 spans the whole domain, but finer
591  // levels do not have such a restriction.
592  // For now, use the bounding box of the boxArray
593  // for normalization, consistent with avg routine.
594 
595  // Bounded box of CC data used for normalization
596  Box bnd_bx = (m_fields[lev][3]->boxArray()).minimalBox();
597 
598  // NOTE: Bounding box must lie on the periodic boundaries
599  // in order to trip the is_per flag
600 
601  // Num components, plane avg, cells per plane
602  Array<int,AMREX_SPACEDIM> is_per = {0,0,0};
603  for (int idim(0); idim < AMREX_SPACEDIM; ++idim) {
604  if ( m_geom[lev].isPeriodic(idim) &&
605  bnd_bx.bigEnd(idim)==domain.bigEnd(idim) &&
606  bnd_bx.smallEnd(idim)==domain.smallEnd(idim) ) { is_per[idim] = 1; }
607  }
608 
609  m_ncell_plane[lev].resize(m_navg);
610  m_plane_average[lev].resize(m_navg);
611  for (int iavg(0); iavg < m_navg; ++iavg) {
612  // Convert bnd_bx to current index type
613  IndexType ixt = m_averages[lev][iavg]->boxArray().ixType();
614  bnd_bx.convert(ixt);
615  IntVect bnd_bx_lo(bnd_bx.loVect());
616  IntVect bnd_bx_hi(bnd_bx.hiVect());
617 
618  m_plane_average[lev][iavg] = zero;
619 
620  m_ncell_plane[lev][iavg] = 1;
621  for (int idim(0); idim < AMREX_SPACEDIM; ++idim) {
622  if (idim != m_face.coordDir()) {
623  if (ixt.nodeCentered(idim) && is_per[idim]) {
624  m_ncell_plane[lev][iavg] *= (bnd_bx_hi[idim] - bnd_bx_lo[idim]);
625  } else {
626  m_ncell_plane[lev][iavg] *= (bnd_bx_hi[idim] - bnd_bx_lo[idim] + 1);
627  }
628  }
629  } // idim
630  } // iavg
631 }

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

Compute number of cells in the region average.

Parameters
[in]levlevel index
132  {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)

Update the rotated fields.

Parameters
[in]levlevel index

Function to set the rotated velocities.

Parameters
[in]levCurrent level
514 {
515  // Peel back the level
516  auto& fields = m_fields[lev];
517  auto& rot_fields = m_rot_fields[lev];
518  auto z_phys_nd = m_z_phys_nd[lev];
519 
520  // Inverse grid size
521  const auto dxInv = m_geom[lev].InvCellSizeArray();
522 
523  // Single MFIter over CC data
524  int imf_cc = 3;
525 
526  // NOTE: The region average reads the rotated velocities out into the ghost
527  // region (see get_ng_fill) and a FillBoundary here would not help an
528  // isolated fine patch, so we rotate the ghost cells as well. We stop
529  // one cell shy of the full ghost region since the rotation reaches one
530  // cell ahead in w (and one node ahead in z_phys).
531  IntVect ngu = max(rot_fields[0]->nGrowVect() - IntVect(1,1,1), IntVect(0));
532  IntVect ngv = max(rot_fields[1]->nGrowVect() - IntVect(1,1,1), IntVect(0));
533 
534  // Populate rotated U & V for terrain
535 #ifdef _OPENMP
536 #pragma omp parallel if (Gpu::notInLaunchRegion())
537 #endif
538  for (MFIter mfi(*fields[imf_cc], TileNoZ()); mfi.isValid(); ++mfi) {
539  Box ubx = mfi.tilebox(IntVect(1,0,0),ngu);
540  Box vbx = mfi.tilebox(IntVect(0,1,0),ngv);
541 
542  const Array4<const Real>& z_phys_arr = z_phys_nd->const_array(mfi);
543 
544  const Array4<const Real>& u_arr = fields[0]->const_array(mfi);
545  const Array4<const Real>& v_arr = fields[1]->const_array(mfi);
546  const Array4<const Real>& w_arr = fields[2]->const_array(mfi);
547 
548  const Array4<Real>& u_rot_arr = rot_fields[0]->array(mfi);
549  const Array4<Real>& v_rot_arr = rot_fields[1]->array(mfi);
550 
551  // U rotated magnitude
552  ParallelFor(ubx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
553  {
554  // Elements of first tangent vector
555  Real met_h_xi = Compute_h_xi_AtIface(i,j,k,dxInv,z_phys_arr);
556  u_rot_arr(i,j,k) = (u_arr(i,j,k) + met_h_xi*w_arr(i,j,k))
557  / std::sqrt(met_h_xi*met_h_xi + one);
558  });
559 
560  // V rotated magnitude
561  ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
562  {
563  // Elements of second tangent vector
564  Real met_h_eta = Compute_h_eta_AtJface(i,j,k,dxInv,z_phys_arr);
565  v_rot_arr(i,j,k) = (v_arr(i,j,k) + met_h_eta*w_arr(i,j,k))
566  / std::sqrt(met_h_eta*met_h_eta + one);
567  });
568  }
569 
570  // Direct copy of other scalar variables
571  MultiFab::Copy(*rot_fields[3],*fields[3],0,0,1,rot_fields[3]->nGrowVect());
572  if (fields[4]) MultiFab::Copy(*rot_fields[4],*fields[4],0,0,1,rot_fields[4]->nGrowVect());
573 }
Compute_h_eta_AtJface
AMREX_GPU_HOST_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:396
Compute_h_xi_AtIface
AMREX_GPU_HOST_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:303

Referenced by compute_averages().

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

◆ set_time_avg_initialized()

void MOSTAverage::set_time_avg_initialized ( const int &  lev)
inline

Declare the time filter at this level to hold meaningful history. This is used on restart, once the filter state has been read back in.

Parameters
[in]levlevel index
290  {
291  if (m_t_avg) { m_t_init[lev] = 1; }
292  }

Referenced by SurfaceLayer::set_mac_avg_initialized().

Here is the caller graph for this function:

◆ set_z_positions_EB()

void MOSTAverage::set_z_positions_EB ( const int &  lev)

Populate embedded-boundary z positions.

Parameters
[in]levlevel index

Function to set K indices for EB.

Parameters
[in]levCurrent level
980 {
981  // See zref_sentinel in ERF_MOSTAverage.H: probe with the sentinel and test the
982  // value, since queryAdd's return value stops meaning "user specified" after the
983  // first parse of the key and this routine runs once per level.
984  Real zref_tmp = zref_sentinel;
985  ParmParse pp(m_pp_prefix);
986  pp.queryAdd("most.zref",zref_tmp);
987  bool read_z = (zref_tmp > Real(0));
988  if (!read_z) { zref_tmp = zref_default; }
989 
990  if (read_z) {
991  m_zref[lev]->setVal( zref_tmp );
992  // Default behavior is to use the first cell center
993  } else {
994  Real m_dz = m_geom[0].CellSize(2);
995  zref_tmp = myhalf * m_dz;
996  m_zref[lev]->setVal( zref_tmp );
997  Print() << "Reference height for MOST set to " << zref_tmp << std::endl;
998  }
999 }

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)

Populate terrain-aware z positions.

Parameters
[in]levlevel index

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

Parameters
[in]levCurrent level
1276 {
1277  AMREX_ALWAYS_ASSERT(m_face.coordDir() == 2);
1278 
1279  // Peel back the level
1280  auto& fields = m_fields[lev];
1281 
1282  // MFIter over CC data
1283  int imf_cc = 3;
1284 
1285  ParmParse pp(m_pp_prefix);
1286  // See zref_sentinel in ERF_MOSTAverage.H: probe with the sentinel and test the
1287  // value, since queryAdd's return value stops meaning "user specified" after the
1288  // first parse of the key and this routine runs once per level.
1289  Real zref_tmp = zref_sentinel;
1290  pp.queryAdd("most.zref",zref_tmp);
1291  bool read_zref = (zref_tmp > Real(0));
1292  if (!read_zref) {
1293  zref_tmp = zref_default;
1294  Print() << "most.zref not specified, query distance default is " << zref_tmp << std::endl;
1295  }
1296  m_zref[lev]->setVal(zref_tmp);
1297  int klo = m_geom[lev].Domain().smallEnd(2);
1298 
1299  // Capture for device
1300  Real d_zref = zref_tmp;
1301  const auto plo = m_geom[lev].ProbLoArray();
1302 
1303  RealVect base;
1304  const auto dx = m_geom[lev].CellSizeArray();
1305  IntVect ng = m_x_pos[lev]->nGrowVect(); ng[2]=0;
1306  const int position_ng = (m_radius > 1) ? m_radius : 1;
1307 
1308  // The positions are read over the averaging stencil at every cell where an
1309  // average is computed, ghost cells included; the box is made nodal so that
1310  // we also cover the U & V averages, which are face centered
1311  const IntVect ng_pos = max(get_ng_fill(lev) + IntVect(m_radius,m_radius,0),
1312  IntVect(position_ng,position_ng,0));
1313 
1314  for (MFIter mfi(*fields[imf_cc], TileNoZ()); mfi.isValid(); ++mfi) {
1315  Box npbx = mfi.tilebox(IntVect(1,1,0),ng_pos);
1316  const Box vbx = mfi.validbox();
1317  Box gtbx = mfi.growntilebox(ng);
1318 
1319  if (vbx.smallEnd(2) != klo) { continue; }
1320 
1321  npbx.makeSlab(2,klo);
1322 
1323  RealBox grb{gtbx,dx.data(),base.dataPtr()};
1324 
1325  const auto z_phys_arr = m_z_phys_nd[lev]->const_array(mfi);
1326  auto x_pos_arr = m_x_pos[lev]->array(mfi);
1327  auto y_pos_arr = m_y_pos[lev]->array(mfi);
1328  auto z_pos_arr = m_z_pos[lev]->array(mfi);
1329  ParallelFor(npbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
1330  {
1331  // Final position at end of vector
1332  x_pos_arr(i,j,k) = plo[0] + ((Real) i + myhalf) * dx[0];
1333  y_pos_arr(i,j,k) = plo[1] + ((Real) j + myhalf) * dx[1];
1334  Real z_bot_face = fourth * ( z_phys_arr(i ,j ,k) + z_phys_arr(i+1,j ,k)
1335  + z_phys_arr(i ,j+1,k) + z_phys_arr(i+1,j+1,k) );
1336  z_pos_arr(i,j,k) = z_bot_face + d_zref;
1337 
1338  // Destination position must be contained on the current process!
1339  Real pos[] = {x_pos_arr(i,j,k)-plo[0],y_pos_arr(i,j,k)-plo[1],myhalf*dx[2]};
1340  amrex::ignore_unused(pos);
1341  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(grb.contains(&pos[0]),
1342  "Query point outside of proc domain!");
1343  });
1344  }
1345 }

Referenced by make_MOSTAverage_at_level().

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

◆ time_avg_is_initialized()

bool MOSTAverage::time_avg_is_initialized ( const int &  lev) const
inline

Return whether the time filter at this level holds meaningful history. Until it does, the next average is taken to be the instantaneous one.

Parameters
[in]levlevel index
279  {
280  return ( m_t_avg && (lev < static_cast<int>(m_t_init.size())) && (m_t_init[lev] == 1) );
281  }

Referenced by SurfaceLayer::mac_avg_is_initialized().

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
[in]zpZ-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
424  {
425  // Search to get z/k
426  bool found = false;
427  int kmax = ubound(z_arr).z;
428  amrex::Real zval = zero;
429  amrex::Real z_target = zp;
430 
431  // Map position to i,j (must be same mapping in cpp file)
432  amrex::Real ireal = (xp - plo[0]) * dxi[0];
433  amrex::Real jreal = (yp - plo[1]) * dxi[1];
434  int i_new = (int) (ireal - myhalf);
435  int j_new = (int) (jreal - myhalf);
436 
437  for (int lk(0); lk<kmax; ++lk) {
438  amrex::Real z_lo = fourth * ( z_arr(i_new,j_new ,lk ) + z_arr(i_new+1,j_new ,lk )
439  + z_arr(i_new,j_new+1,lk ) + z_arr(i_new+1,j_new+1,lk ) );
440  amrex::Real z_hi = fourth * ( z_arr(i_new,j_new ,lk+1) + z_arr(i_new+1,j_new ,lk+1)
441  + z_arr(i_new,j_new+1,lk+1) + z_arr(i_new+1,j_new+1,lk+1) );
442  if (in_cell_z(z_target, z_lo, z_hi)) {
443  found = true;
444  // Weight between the two cell centers that bracket z_target by
445  // their physical heights, so the interpolation is linear in z
446  // on a stretched or terrain-following column. A weight taken
447  // from the fraction of the containing cell alone (as before)
448  // is exact only for equal cell heights: at the 10 m face of a
449  // 10 m / 11 m pair it returned the wind at 10.25 m. A center
450  // beyond the nodes searched here is mirrored about the shared
451  // face, which gives the equal-height weight at those ends.
452  amrex::Real zc = myhalf * (z_lo + z_hi);
453  if (z_target >= zc) {
454  amrex::Real zc_up = two * z_hi - zc;
455  if (lk + 2 <= kmax) {
456  amrex::Real z_up = fourth * ( z_arr(i_new,j_new ,lk+2) + z_arr(i_new+1,j_new ,lk+2)
457  + z_arr(i_new,j_new+1,lk+2) + z_arr(i_new+1,j_new+1,lk+2) );
458  zc_up = myhalf * (z_hi + z_up);
459  }
460  zval = (amrex::Real) (lk + 1) + (z_target - zc) / (zc_up - zc);
461  } else {
462  amrex::Real zc_dn = two * z_lo - zc;
463  if (lk >= 1) {
464  amrex::Real z_dn = fourth * ( z_arr(i_new,j_new ,lk-1) + z_arr(i_new+1,j_new ,lk-1)
465  + z_arr(i_new,j_new+1,lk-1) + z_arr(i_new+1,j_new+1,lk-1) );
466  zc_dn = myhalf * (z_dn + z_lo);
467  }
468  zval = (amrex::Real) lk + (z_target - zc_dn) / (zc - zc_dn);
469  }
470  break;
471  }
472  }
473 
474  amrex::ignore_unused(found);
475  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(found, "MOSTAverage: Height above terrain not found, try increasing z_ref!");
476 
477  // NOTE: This is the point ahead of the current i,j (e.g. i/j_new + 1);
478  // in z, floor(zval) is the upper of the two bracketing cells
479  const amrex::RealVect lx(ireal + myhalf, jreal + myhalf, zval);
480 
481  const amrex::IntVect ijk = lx.floor();
482 
483  int i = ijk[0]; int j = ijk[1]; int k = ijk[2];
484 
485  // Convert ijk (IntVect) to a RealVect explicitly
486  amrex::RealVect ijk_r(static_cast<amrex::Real>(ijk[0]),
487  static_cast<amrex::Real>(ijk[1]),
488  static_cast<amrex::Real>(ijk[2]));
489 
490  // Weights
491  const amrex::RealVect sx_hi = lx - ijk_r;
492  const amrex::RealVect sx_lo = one - sx_hi;
493 
494  for (int n = 0; n < interp_comp; n++) {
495  interp_vals[n] = sx_lo[0]*sx_lo[1]*sx_lo[2]*interp_array(i-1, j-1, k-1,n) +
496  sx_lo[0]*sx_lo[1]*sx_hi[2]*interp_array(i-1, j-1, k ,n) +
497  sx_lo[0]*sx_hi[1]*sx_lo[2]*interp_array(i-1, j , k-1,n) +
498  sx_lo[0]*sx_hi[1]*sx_hi[2]*interp_array(i-1, j , k ,n) +
499  sx_hi[0]*sx_lo[1]*sx_lo[2]*interp_array(i , j-1, k-1,n) +
500  sx_hi[0]*sx_lo[1]*sx_hi[2]*interp_array(i , j-1, k ,n) +
501  sx_hi[0]*sx_hi[1]*sx_lo[2]*interp_array(i , j , k-1,n) +
502  sx_hi[0]*sx_hi[1]*sx_hi[2]*interp_array(i , j , k ,n);
503  }
504  }
lx
const Real lx
Definition: ERF_InitCustomPertVels_CloudChamber.H:10
zc
const Real zc
Definition: ERF_InitCustomPert_ABL.H:56
two
constexpr amrex::Real two
Definition: ERF_NumericalConstants.H:38

Referenced by compute_plane_averages(), and compute_region_averages().

Here is the call graph for this function:
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 
)

Reset pointers to field MultiFabs.

Parameters
[in]levlevel index
[in]vars_oldold-time state and velocity fields
[in]Theta_primprimitive potential-temperature fields by level
[in]Qv_primprimitive water-vapor fields by level
[in]Qr_primprimitive rain-water fields by level

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
[in]Qv_primPrimitive water-vapor component at each level
[in]Qr_primPrimitive rain-water component at each level
447 {
448  m_fields[lev][0] = &vars_old[lev][Vars::xvel];
449  m_fields[lev][1] = &vars_old[lev][Vars::yvel];
450  m_fields[lev][2] = &vars_old[lev][Vars::zvel];
451  m_fields[lev][3] = Theta_prim[lev].get();
452  m_fields[lev][4] = Qv_prim[lev].get();
453  m_fields[lev][5] = Qr_prim[lev].get();
454 }

Referenced by SurfaceLayer::update_mac_ptrs().

Here is the caller graph for this function:

◆ validate_lateral_reference_stencil()

void MOSTAverage::validate_lateral_reference_stencil ( const int &  lev)

Verify that the lateral face-owned FAB contains the reference stencil.

Parameters
[in]levlevel index

This is a conservative preflight check for the direct regional-average path. It verifies that the reference point plus averaging radius, with allowance for staggered-field neighbors, is available in the owned FAB.

Verify that lateral face-owned FABs contain the direct regional stencil.

Parameters
[in]levCurrent level.

The check uses the full cell-centered source field as a conservative layout proxy for the staggered fields used by the regional kernels. It includes ghost cells and one extra cell on each high side for staggered neighbors.

TODO: Remove this check if a distributed reference-plane staging implementation is added to the direct lateral regional path.

941 {
942  const int dir = m_face.coordDir();
943  if (dir == 2 || m_interp || m_norm_vec || !m_k_indx[lev]) { return; }
944 
945  const MultiFab& source = *m_fields[lev][3];
946  const int face_index = m_face.isLow()
947  ? m_geom[lev].Domain().smallEnd(dir)
948  : m_geom[lev].Domain().bigEnd(dir);
949  const int reference = m_k_indx[lev]->min(0);
950 
951  for (MFIter mfi(source, false); mfi.isValid(); ++mfi) {
952  const Box& valid_box = mfi.validbox();
953  const bool owns_face = m_face.isLow()
954  ? valid_box.smallEnd(dir) == face_index
955  : valid_box.bigEnd(dir) == face_index;
956  if (!owns_face) { continue; }
957 
958  Box required = valid_box;
959  required.setSmall(dir, reference - m_radius);
960  required.setBig(dir, reference + m_radius);
961  required.growHi(0, 1);
962  required.growHi(1, 1);
963  required.growHi(2, 1);
964 
965  Box available = valid_box;
966  available.grow(source.nGrowVect());
967  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
968  available.contains(required),
969  "MOST lateral reference plus averaging radius is outside the face-owned FAB.");
970  }
971 }

Referenced by make_MOSTAverage_at_level().

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

◆ write_averages()

void MOSTAverage::write_averages ( const int &  lev)

Write averages on the 2D MultiFab.

Parameters
[in]levlevel index

Function to write averages to text file.

Parameters
[in]levCurrent level
3222 {
3223  // Peel back the level
3224  auto& fields = m_fields[lev];
3225  auto& averages = m_averages[lev];
3226 
3227  // MFIter on CC
3228  int imf_cc = 3;
3229 
3230 
3231  int navg = m_navg - 1;
3232 
3233  std::ofstream ofile;
3234  ofile.open (std::string("MOST_averages_" + std::to_string(m_face) + ".txt"));
3235  ofile << "Averages computed via MOSTAverages class:\n";
3236 
3237  const int dir = m_face.coordDir();
3238  int sm_index;
3239  if (m_face.isLow()) {
3240  sm_index = m_geom[lev].Domain().smallEnd(dir);
3241  } else {
3242  sm_index = m_geom[lev].Domain().bigEnd(dir);
3243  }
3244 
3245  for (MFIter mfi(*fields[imf_cc], TileNoZ()); mfi.isValid(); ++mfi) {
3246  Box pbx = mfi.tilebox();
3247 
3248  if (m_face.isLow()) {
3249  if (mfi.validbox().smallEnd(dir) != sm_index ||
3250  pbx.smallEnd(dir) != sm_index) {
3251  continue;
3252  }
3253  } else {
3254  if (mfi.validbox().bigEnd(dir) != sm_index ||
3255  pbx.bigEnd(dir) != sm_index) {
3256  continue;
3257  }
3258  }
3259 
3260  pbx.setSmall(dir, sm_index); pbx.setBig(dir, sm_index);
3261 
3262  int il = pbx.smallEnd(0); int iu = pbx.bigEnd(0);
3263  int jl = pbx.smallEnd(1); int ju = pbx.bigEnd(1);
3264  int kl = pbx.smallEnd(2); int ku = pbx.bigEnd(2);
3265 
3266  for (int k(kl); k <= ku; ++k) {
3267  for (int j(jl); j <= ju; ++j) {
3268  for (int i(il); i <= iu; ++i) {
3269  if (dir == 0) {
3270  ofile << "(J,K): " << "(" << j << "," << k << ")" << "\n";
3271  } else if (dir == 1) {
3272  ofile << "(I,K): " << "(" << i << "," << k << ")" << "\n";
3273  } else {
3274  ofile << "(I,J): " << "(" << i << "," << j << ")" << "\n";
3275  }
3276 
3277  for (int iavg(0); iavg <= navg; ++iavg) {
3278  auto mf_arr = averages[iavg]->array(mfi);
3279  ofile << "iavg val: "
3280  << iavg << ' '
3281  << mf_arr(i,j,k) << "\n";
3282  }
3283  ofile << "\n";
3284  }
3285  }
3286  }
3287  }
3288  ofile.close();
3289 }
Here is the call graph for this function:

◆ write_k_indices()

void MOSTAverage::write_k_indices ( const int &  lev)

Write k-index data.

Parameters
[in]levlevel index

Function to write the K indices to text file.

Parameters
[in]levCurrent level
3024 {
3025  // Peel back the level
3026  auto& fields = m_fields[lev];
3027  auto& k_indx = m_k_indx[lev];
3028 
3029  // MFIter on CC
3030  int imf_cc = 3;
3031 
3032  std::ofstream ofile;
3033  ofile.open ("MOST_k_indices_" + std::to_string(m_face) + ".txt");
3034  ofile << "K indices used to compute averages via MOSTAverages class:\n";
3035 
3036  const int dir = m_face.coordDir();
3037  int sm_index;
3038  if (m_face.isLow()) {
3039  sm_index = m_geom[lev].Domain().smallEnd(dir);
3040  } else {
3041  sm_index = m_geom[lev].Domain().bigEnd(dir);
3042  }
3043 
3044  for (MFIter mfi(*fields[imf_cc], TileNoZ()); mfi.isValid(); ++mfi) {
3045  Box pbx = mfi.tilebox();
3046 
3047  if (m_face.isLow()) {
3048  if (mfi.validbox().smallEnd(dir) != sm_index ||
3049  pbx.smallEnd(dir) != sm_index) {
3050  continue;
3051  }
3052  } else {
3053  if (mfi.validbox().bigEnd(dir) != sm_index ||
3054  pbx.bigEnd(dir) != sm_index) {
3055  continue;
3056  }
3057  }
3058 
3059  pbx.setSmall(dir, sm_index); pbx.setBig(dir, sm_index);
3060 
3061  int il = pbx.smallEnd(0); int iu = pbx.bigEnd(0);
3062  int jl = pbx.smallEnd(1); int ju = pbx.bigEnd(1);
3063  int kl = pbx.smallEnd(2); int ku = pbx.bigEnd(2);
3064 
3065  auto k_arr = k_indx->array(mfi);
3066  for (int k(kl); k <= ku; ++k) {
3067  for (int j(jl); j <= ju; ++j) {
3068  for (int i(il); i <= iu; ++i) {
3069  if (dir == 0) {
3070  ofile << "(J,K): " << "(" << j << "," << k << ")" << "\n";
3071  } else if (dir == 1) {
3072  ofile << "(I,K): " << "(" << i << "," << k << ")" << "\n";
3073  } else {
3074  ofile << "(I,J): " << "(" << i << "," << j << ")" << "\n";
3075  }
3076  ofile << "K_ind: "
3077  << k_arr(i,j,k) << "\n";
3078  ofile << "\n";
3079  }
3080  }
3081  }
3082  }
3083  ofile.close();
3084 }
Here is the call graph for this function:

◆ write_norm_indices()

void MOSTAverage::write_norm_indices ( const int &  lev)

Write normal-index data.

Parameters
[in]levlevel index

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

Parameters
[in]levCurrent level
3094 {
3095  // Peel back the level
3096  auto& fields = m_fields[lev];
3097  auto& k_indx = m_k_indx[lev];
3098  auto& j_indx = m_j_indx[lev];
3099  auto& i_indx = m_i_indx[lev];
3100 
3101  // MFIter on CC
3102  int imf_cc = 3;
3103 
3104  std::ofstream ofile;
3105  ofile.open ("MOST_ijk_indices_" + std::to_string(m_face) + ".txt");
3106  ofile << "IJK indices used to compute averages via MOSTAverages class:\n";
3107 
3108  const int dir = m_face.coordDir();
3109  int sm_index;
3110  if (m_face.isLow()) {
3111  sm_index = m_geom[lev].Domain().smallEnd(dir);
3112  } else {
3113  sm_index = m_geom[lev].Domain().bigEnd(dir);
3114  }
3115 
3116  for (MFIter mfi(*fields[imf_cc], TileNoZ()); mfi.isValid(); ++mfi) {
3117  Box pbx = mfi.tilebox();
3118 
3119  if (m_face.isLow()) {
3120  if (mfi.validbox().smallEnd(dir) != sm_index ||
3121  pbx.smallEnd(dir) != sm_index) {
3122  continue;
3123  }
3124  } else {
3125  if (mfi.validbox().bigEnd(dir) != sm_index ||
3126  pbx.bigEnd(dir) != sm_index) {
3127  continue;
3128  }
3129  }
3130 
3131  pbx.setSmall(dir, sm_index); pbx.setBig(dir, sm_index);
3132 
3133  int il = pbx.smallEnd(0); int iu = pbx.bigEnd(0);
3134  int jl = pbx.smallEnd(1); int ju = pbx.bigEnd(1);
3135  int kl = pbx.smallEnd(2); int ku = pbx.bigEnd(2);
3136 
3137  auto k_arr = k_indx->array(mfi);
3138  auto j_arr = j_indx ? j_indx->array(mfi) : Array4<int> {};
3139  auto i_arr = i_indx ? i_indx->array(mfi) : Array4<int> {};
3140 
3141  for (int k(kl); k <= ku; ++k) {
3142  for (int j(jl); j <= ju; ++j) {
3143  for (int i(il); i <= iu; ++i) {
3144  ofile << "(I1,J1,K1): " << "(" << i << "," << j << "," << k << ")" << "\n";
3145 
3146  // This diagnostic reports the cell-centered reference used
3147  // for the CC fields. The high-face +1 adjustment applies
3148  // only to the face-centered velocity normal to the wall.
3149  const int ref = k_arr(i,j,k);
3150  int im = i_arr ? i_arr(i,j,k) : i;
3151  int jm = j_arr ? j_arr(i,j,k) : j;
3152  int km = k;
3153  if (dir == 0) {
3154  im = ref;
3155  } else if (dir == 1) {
3156  jm = ref;
3157  } else {
3158  km = ref;
3159  }
3160 
3161  ofile << "(I2,J2,K2): "
3162  << "(" << im << "," << jm << "," << km << ")" << "\n";
3163  ofile << "\n";
3164  }
3165  }
3166  }
3167  }
3168  ofile.close();
3169 }
Here is the call graph for this function:

◆ write_xz_positions()

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

Write XZ planar positions.

Parameters
[in]levlevel index
[in]jy-index of the plane to write

Function to write X & Z positions to text file.

Parameters
[in]levCurrent level
[in]jIndex in y-dir
3181 {
3182  // Peel back the level
3183  auto& fields = m_fields[lev];
3184  auto& x_pos_mf = m_x_pos[lev];
3185  auto& z_pos_mf = m_z_pos[lev];
3186 
3187  // MFIter on CC
3188  int imf_cc = 3;
3189 
3190  int klo = m_geom[lev].Domain().smallEnd(2);
3191 
3192  std::ofstream ofile;
3193  ofile.open ("MOST_xz_positions.txt");
3194 
3195  for (MFIter mfi(*fields[imf_cc], TileNoZ()); mfi.isValid(); ++mfi) {
3196  Box pbx = mfi.tilebox();
3197 
3198  if(pbx.smallEnd(2) != klo) { continue; }
3199 
3200  pbx.makeSlab(2,klo);
3201 
3202  int il = pbx.smallEnd(0); int iu = pbx.bigEnd(0);
3203 
3204  auto x_pos_arr = x_pos_mf->array(mfi);
3205  auto z_pos_arr = z_pos_mf->array(mfi);
3206 
3207  int k = 0;
3208  for (int i(il); i <= iu; ++i)
3209  ofile << x_pos_arr(i,j,k) << ' ' << z_pos_arr(i,j,k) << "\n";
3210  }
3211  ofile.close();
3212 }
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(), extrap_ghost_cells(), get_average(), get_ng_fill(), 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_face

amrex::Orientation MOSTAverage::m_face
protected

Referenced by compute_eb_averages(), compute_plane_averages(), compute_region_averages(), extrap_ghost_cells(), fill_planar_boundary(), get_ng_fill(), make_MOSTAverage_at_level(), set_k_indices_N(), set_k_indices_T(), set_norm_indices_T(), set_norm_positions_T(), set_plane_normalization(), set_z_positions_T(), validate_lateral_reference_stencil(), write_averages(), write_k_indices(), and write_norm_indices().

◆ 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(), get_ng_fill(), 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(), validate_lateral_reference_stencil(), 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(), fill_planar_boundary(), 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(), validate_lateral_reference_stencil(), 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(), make_MOSTAverage_at_level(), and validate_lateral_reference_stencil().

◆ 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(), validate_lateral_reference_stencil(), 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

Referenced by set_k_indices_N().

◆ m_navg

int MOSTAverage::m_navg {9}
protected

Referenced by compute_eb_averages(), compute_plane_averages(), compute_region_averages(), get_navg(), get_ng_fill(), 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 compute_plane_averages(), compute_region_averages(), make_MOSTAverage_at_level(), and validate_lateral_reference_stencil().

◆ m_nvar

int MOSTAverage::m_nvar {6}
protected

Referenced by get_ng_fill(), and make_MOSTAverage_at_level().

◆ m_planar_bndry

amrex::Vector<PlanarBoundary> MOSTAverage::m_planar_bndry
protected

Referenced by fill_planar_boundary(), and 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(), get_plane_average(), set_eb_normalization(), set_plane_average(), 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(), get_ng_fill(), set_k_indices_N(), set_k_indices_T(), set_norm_indices_T(), set_norm_positions_T(), set_region_normalization(), set_z_positions_T(), and validate_lateral_reference_stencil().

◆ 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(), do_time_averaging(), make_MOSTAverage_at_level(), set_time_avg_initialized(), and time_avg_is_initialized().

◆ m_t_init

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

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

◆ m_terrain_type

TerrainType MOSTAverage::m_terrain_type
protected

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

◆ m_time_window

double MOSTAverage::m_time_window {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_zlevels_stag

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

Referenced by set_k_indices_N().

◆ 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().

◆ zref_sentinel

const amrex::Real MOSTAverage::zref_sentinel = amrex::Real(-1.0)
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().

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