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)
 
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:34

Referenced by k_index_below(), and 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
1429 {
1430  if (m_rotate) set_rotated_fields(lev);
1431 
1432  switch(m_policy) {
1433  case 0: // Standard plane average
1434  compute_plane_averages(lev);
1435  break;
1436  case 1: // Local region/point
1437  compute_region_averages(lev);
1438  break;
1439  case 2: // EB Terrain average
1440  compute_eb_averages(lev);
1441  break;
1442  default:
1443  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(false, "Unknown policy for MOSTAverage!");
1444  }
1445 
1446  // We have initialized the averages
1447  if (m_t_avg) m_t_init[lev] = 1;
1448 }
AMREX_ALWAYS_ASSERT_WITH_MESSAGE
AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_cloud_chamber_config.active, "Cloud Chamber: initializer reached without a parsed configuration")
MOSTAverage::m_t_avg
bool m_t_avg
Definition: ERF_MOSTAverage.H:593
MOSTAverage::compute_plane_averages
void compute_plane_averages(const int &lev)
Definition: ERF_MOSTAverage.cpp:1457
MOSTAverage::m_policy
int m_policy
Definition: ERF_MOSTAverage.H:557
MOSTAverage::compute_region_averages
void compute_region_averages(const int &lev)
Definition: ERF_MOSTAverage.cpp:2056
MOSTAverage::m_t_init
amrex::Vector< int > m_t_init
Definition: ERF_MOSTAverage.H:594
MOSTAverage::set_rotated_fields
void set_rotated_fields(const int &lev)
Definition: ERF_MOSTAverage.cpp:512
MOSTAverage::m_rotate
bool m_rotate
Definition: ERF_MOSTAverage.H:558
MOSTAverage::compute_eb_averages
void compute_eb_averages(const int &lev)
Definition: ERF_MOSTAverage.cpp:2689
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
2690 {
2691  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
2692  m_face.coordDir() == 2,
2693  "EB MOST averaging currently supports only z faces.");
2694  AMREX_ALWAYS_ASSERT(m_eb_vec[lev] != nullptr);
2695 
2696  // Peel back the level
2697  auto& fields = m_fields[lev];
2698  auto& averages = m_averages[lev];
2699  auto& plane_average = m_plane_average[lev];
2700 
2701  // Get EB data
2702  const auto& cc_flags = m_eb_vec[lev]->get_const_factory()->getMultiEBCellFlagFab();
2703  auto cc_afrac = m_eb_vec[lev]->get_const_factory()->getAreaFrac();
2704  const auto& cc_bnorm = m_eb_vec[lev]->get_const_factory()->getBndryNormal();
2705  const auto& u_vfrac = m_eb_vec[lev]->get_u_const_factory()->getVolFrac();
2706  const auto& v_vfrac = m_eb_vec[lev]->get_v_const_factory()->getVolFrac();
2707  const auto& w_vfrac = m_eb_vec[lev]->get_w_const_factory()->getVolFrac();
2708 
2709  // Get geometry for cell sizes
2710  auto const& dx_arr = m_geom[lev].CellSizeArray();
2711  Real dx = dx_arr[0];
2712  Real dy = dx_arr[1];
2713  Real dz = dx_arr[2];
2714 
2715  // Set factors for time averaging
2716  Real d_fact_new, d_fact_old;
2717  if (m_t_avg && m_t_init[lev]) {
2718  d_fact_new = m_fact_new;
2719  d_fact_old = m_fact_old;
2720  } else {
2721  d_fact_new = one;
2722  d_fact_old = zero;
2723  }
2724 
2725  // GPU array to accumulate averages into
2726  Gpu::DeviceVector<Real> pavg(plane_average.size(), zero);
2727  Real* plane_avg = pavg.data();
2728 
2729  // Vectors for normalization and buffer storage
2730  Vector<Real> denom(plane_average.size(),zero);
2731  Vector<Real> val_old(plane_average.size(),zero);
2732 
2733  //
2734  //----------------------------------------------------------
2735  // Averages for U, V, and tangential velocity (in local coordinate)
2736  //----------------------------------------------------------
2737  //
2738  {
2739  denom[0] = one / m_total_bndry_area[lev][0];
2740  val_old[0] = plane_average[0]*d_fact_old;
2741  denom[1] = one / m_total_bndry_area[lev][1];
2742  val_old[1] = plane_average[1]*d_fact_old;
2743 
2744  int iavg = m_navg - 1; // Tangential velocity magnitude
2745  denom[iavg] = one / m_total_bndry_area[lev][iavg];
2746  val_old[iavg] = plane_average[iavg]*d_fact_old;
2747 
2748  const Real Vsg = m_Vsg[lev]; // Subgrid scale velocity
2749 
2750  //
2751  // NOTE: deliberately not threaded. The body is a pure reduction into the shared
2752  // plane_avg entries, and on the host Gpu::deviceReduceSum is an unordered
2753  // "#pragma omp atomic" add. Threading it would make the surface-layer averages
2754  // -- and through them the surface fluxes -- depend on thread scheduling in the
2755  // last bits. See the same note in PlaneAverage::compute_averages.
2756  //
2757  for (MFIter mfi(*fields[3], TileNoZ()); mfi.isValid(); ++mfi) {
2758  const auto& flag = cc_flags[mfi];
2759 
2760  // Skip boxes that are not singlevalued (MultiCutFab only has data for singlevalued boxes)
2761  if (flag.getType() != FabType::singlevalued) continue;
2762 
2763  Box bx = mfi.tilebox(); // Full 3D box
2764 
2765  // Get EB arrays
2766  auto const flag_arr = flag.const_array();
2767  auto const afrac_x = cc_afrac[0]->const_array(mfi);
2768  auto const afrac_y = cc_afrac[1]->const_array(mfi);
2769  auto const afrac_z = cc_afrac[2]->const_array(mfi);
2770  auto const bnorm_arr = cc_bnorm.const_array(mfi);
2771  auto const u_vf_arr = u_vfrac.const_array(mfi);
2772  auto const v_vf_arr = v_vfrac.const_array(mfi);
2773  auto const w_vf_arr = w_vfrac.const_array(mfi);
2774 
2775  // Get velocity arrays
2776  auto const u_arr = fields[0]->const_array(mfi);
2777  auto const v_arr = fields[1]->const_array(mfi);
2778  auto const w_arr = fields[2]->const_array(mfi);
2779 
2780  ParallelFor(Gpu::KernelInfo().setReduction(true), bx, [=]
2781  AMREX_GPU_DEVICE(int i, int j, int k, Gpu::Handler const& handler) noexcept
2782  {
2783  // Area-weighted averaging over cut cells at any k
2784  if (flag_arr(i,j,k).isSingleValued()) {
2785  // Compute area from face-centered area fractions
2786  Real axm = afrac_x(i ,j ,k );
2787  Real axp = afrac_x(i+1,j ,k );
2788  Real aym = afrac_y(i ,j ,k );
2789  Real ayp = afrac_y(i ,j+1,k );
2790  Real azm = afrac_z(i ,j ,k );
2791  Real azp = afrac_z(i ,j ,k+1);
2792 
2793  Real adx = (axm - axp) * dy * dz;
2794  Real ady = (aym - ayp) * dx * dz;
2795  Real adz = (azm - azp) * dx * dy;
2796 
2797  Real area = std::sqrt(adx*adx + ady*ady + adz*adz);
2798 
2799  // Volume-weighted interpolation of velocities to cell center
2800  Real vf_u_lo = u_vf_arr(i,j,k);
2801  Real vf_u_hi = u_vf_arr(i+1,j,k);
2802  Real sum_vf_u = vf_u_lo + vf_u_hi;
2803  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;
2804 
2805  Real vf_v_lo = v_vf_arr(i,j,k);
2806  Real vf_v_hi = v_vf_arr(i,j+1,k);
2807  Real sum_vf_v = vf_v_lo + vf_v_hi;
2808  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;
2809 
2810  Real vf_w_lo = w_vf_arr(i,j,k);
2811  Real vf_w_hi = w_vf_arr(i,j,k+1);
2812  Real sum_vf_w = vf_w_lo + vf_w_hi;
2813  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;
2814 
2815  // Get normal vector components
2816  Real nx = bnorm_arr(i,j,k,0);
2817  Real ny = bnorm_arr(i,j,k,1);
2818  Real nz = bnorm_arr(i,j,k,2);
2819 
2820  // Compute tangential velocity components
2821  Real v_dot_n = u_cc*nx + v_cc*ny + w_cc*nz;
2822  Real u_tangent = u_cc - v_dot_n * nx;
2823  Real v_tangent = v_cc - v_dot_n * ny;
2824  Real mag = std::sqrt(u_tangent*u_tangent + v_tangent*v_tangent + Vsg*Vsg);
2825 
2826  // Area-weighted sum
2827  Real val_u = u_tangent * area;
2828  Real val_v = v_tangent * area;
2829  Real val_mag = mag * area;
2830 
2831  Gpu::deviceReduceSum(&plane_avg[0], val_u, handler);
2832  Gpu::deviceReduceSum(&plane_avg[1], val_v, handler);
2833  Gpu::deviceReduceSum(&plane_avg[iavg], val_mag, handler);
2834  }
2835  });
2836  }
2837  }
2838 
2839  //
2840  //----------------------------------------------------------
2841  // Averages for T,Qv (cell-centered scalars)
2842  //----------------------------------------------------------
2843  //
2844  for (int imf(3); imf < 5; ++imf) {
2845 
2846  // Continue if no valid Qv pointer
2847  if (!fields[imf]) continue;
2848 
2849  denom[imf] = one / m_total_bndry_area[lev][imf];
2850  val_old[imf] = plane_average[imf]*d_fact_old;
2851 
2852  //
2853  // NOTE: deliberately not threaded. The body is a pure reduction into the shared
2854  // plane_avg entries, and on the host Gpu::deviceReduceSum is an unordered
2855  // "#pragma omp atomic" add. Threading it would make the surface-layer averages
2856  // -- and through them the surface fluxes -- depend on thread scheduling in the
2857  // last bits. See the same note in PlaneAverage::compute_averages.
2858  //
2859  for (MFIter mfi(*fields[imf], TileNoZ()); mfi.isValid(); ++mfi) {
2860  const auto& flag = cc_flags[mfi];
2861 
2862  // Skip boxes that are not singlevalued (MultiCutFab only has data for singlevalued boxes)
2863  if (flag.getType() != FabType::singlevalued) continue;
2864 
2865  Box bx = mfi.tilebox(); // Full 3D box
2866  auto const flag_arr = flag.const_array();
2867  auto const afrac_x = cc_afrac[0]->const_array(mfi);
2868  auto const afrac_y = cc_afrac[1]->const_array(mfi);
2869  auto const afrac_z = cc_afrac[2]->const_array(mfi);
2870  auto const mf_arr = fields[imf]->const_array(mfi);
2871 
2872  ParallelFor(Gpu::KernelInfo().setReduction(true), bx, [=]
2873  AMREX_GPU_DEVICE(int i, int j, int k, Gpu::Handler const& handler) noexcept
2874  {
2875  // Area-weighted averaging over cut cells at any k
2876  if (flag_arr(i,j,k).isSingleValued()) {
2877  // Compute area from face-centered area fractions
2878  Real axm = afrac_x(i ,j ,k );
2879  Real axp = afrac_x(i+1,j ,k );
2880  Real aym = afrac_y(i ,j ,k );
2881  Real ayp = afrac_y(i ,j+1,k );
2882  Real azm = afrac_z(i ,j ,k );
2883  Real azp = afrac_z(i ,j ,k+1);
2884 
2885  Real adx = (axm - axp) * dy * dz;
2886  Real ady = (aym - ayp) * dx * dz;
2887  Real adz = (azm - azp) * dx * dy;
2888 
2889  Real area = std::sqrt(adx*adx + ady*ady + adz*adz);
2890 
2891  Real val = mf_arr(i,j,k) * area;
2892  Gpu::deviceReduceSum(&plane_avg[imf], val, handler);
2893  }
2894  });
2895  }
2896  }
2897 
2898  //
2899  //------------------------------------------------------------------------
2900  // Averages for virtual potential temperature
2901  //------------------------------------------------------------------------
2902  //
2903  if (fields[4]) // We have water vapor
2904  {
2905  int iavg = 5;
2906  denom[iavg] = one / m_total_bndry_area[lev][iavg];
2907  val_old[iavg] = plane_average[iavg]*d_fact_old;
2908 
2909  //
2910  // NOTE: deliberately not threaded. The body is a pure reduction into the shared
2911  // plane_avg entries, and on the host Gpu::deviceReduceSum is an unordered
2912  // "#pragma omp atomic" add. Threading it would make the surface-layer averages
2913  // -- and through them the surface fluxes -- depend on thread scheduling in the
2914  // last bits. See the same note in PlaneAverage::compute_averages.
2915  //
2916  for (MFIter mfi(*fields[4], TileNoZ()); mfi.isValid(); ++mfi)
2917  {
2918  const auto& flag = cc_flags[mfi];
2919 
2920  // Skip boxes that are not singlevalued (MultiCutFab only has data for singlevalued boxes)
2921  if (flag.getType() != FabType::singlevalued) continue;
2922 
2923  Box bx = mfi.tilebox(); // Full 3D box
2924  auto const flag_arr = flag.const_array();
2925  auto const afrac_x = cc_afrac[0]->const_array(mfi);
2926  auto const afrac_y = cc_afrac[1]->const_array(mfi);
2927  auto const afrac_z = cc_afrac[2]->const_array(mfi);
2928 
2929  const Array4<Real const> T_mf_arr = fields[3]->const_array(mfi);
2930  const Array4<Real const> qv_mf_arr = fields[4]->const_array(mfi);
2931  const Array4<Real const> qr_mf_arr = (fields[5]) ? fields[5]->const_array(mfi) :
2932  Array4<const Real> {};
2933 
2934  ParallelFor(Gpu::KernelInfo().setReduction(true), bx, [=]
2935  AMREX_GPU_DEVICE(int i, int j, int k, Gpu::Handler const& handler) noexcept
2936  {
2937  if (flag_arr(i,j,k).isSingleValued()) {
2938  // Compute area from face-centered area fractions
2939  Real axm = afrac_x(i ,j ,k );
2940  Real axp = afrac_x(i+1,j ,k );
2941  Real aym = afrac_y(i ,j ,k );
2942  Real ayp = afrac_y(i ,j+1,k );
2943  Real azm = afrac_z(i ,j ,k );
2944  Real azp = afrac_z(i ,j ,k+1);
2945 
2946  Real adx = (axm - axp) * dy * dz;
2947  Real ady = (aym - ayp) * dx * dz;
2948  Real adz = (azm - azp) * dx * dy;
2949 
2950  Real area = std::sqrt(adx*adx + ady*ady + adz*adz);
2951 
2952  Real vfac;
2953  if (qr_mf_arr) {
2954  // We also have liquid water
2955  vfac = one + epsv*qv_mf_arr(i,j,k) - qr_mf_arr(i,j,k);
2956  } else {
2957  vfac = one + epsv*qv_mf_arr(i,j,k);
2958  }
2959  const Real val = T_mf_arr(i,j,k) * vfac * area;
2960  Gpu::deviceReduceSum(&plane_avg[iavg], val, handler);
2961  }
2962  });
2963  }
2964  }
2965  else // copy temperature
2966  {
2967  int iavg = m_navg - 4;
2968  denom[iavg] = one / m_total_bndry_area[lev][iavg];
2969  // plane_avg[iavg] = plane_avg[2]
2970  Gpu::copy(Gpu::deviceToDevice, pavg.begin() + 3, pavg.begin() + 4,
2971  pavg.begin() + iavg);
2972  }
2973 
2974  // Copy to host and sum across procs
2975  Gpu::copy(Gpu::deviceToHost, pavg.begin(), pavg.end(), plane_average.begin());
2976  ParallelDescriptor::ReduceRealSum(plane_average.data(), plane_average.size());
2977 
2978  // Normalize by total area and apply time averaging
2979  for (int iavg(0); iavg < m_navg; ++iavg){
2980  plane_average[iavg] *= denom[iavg]*d_fact_new;
2981  plane_average[iavg] += val_old[iavg];
2982  averages[iavg]->setVal(plane_average[iavg]);
2983  }
2984 }
epsv
constexpr amrex::Real epsv
Definition: ERF_Constants.H:40
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:30
zero
constexpr amrex::Real zero
Definition: ERF_NumericalConstants.H:29
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:555
MOSTAverage::m_averages
amrex::Vector< amrex::Vector< std::unique_ptr< amrex::MultiFab > > > m_averages
Definition: ERF_MOSTAverage.H:566
MOSTAverage::m_total_bndry_area
amrex::Vector< amrex::Vector< amrex::Real > > m_total_bndry_area
Definition: ERF_MOSTAverage.H:589
MOSTAverage::m_face
amrex::Orientation m_face
Definition: ERF_MOSTAverage.H:543
MOSTAverage::m_Vsg
amrex::Vector< amrex::Real > m_Vsg
Definition: ERF_MOSTAverage.H:601
MOSTAverage::m_plane_average
amrex::Vector< amrex::Vector< amrex::Real > > m_plane_average
Definition: ERF_MOSTAverage.H:573
MOSTAverage::m_fact_new
amrex::Real m_fact_new
Definition: ERF_MOSTAverage.H:596
MOSTAverage::m_fields
amrex::Vector< amrex::Vector< amrex::MultiFab * > > m_fields
Definition: ERF_MOSTAverage.H:545
MOSTAverage::m_eb_vec
amrex::Vector< const eb_ * > m_eb_vec
Definition: ERF_MOSTAverage.H:588
MOSTAverage::m_fact_old
amrex::Real m_fact_old
Definition: ERF_MOSTAverage.H:596
MOSTAverage::m_geom
const amrex::Vector< amrex::Geometry > m_geom
Definition: ERF_MOSTAverage.H:544
WDM6SedCellScratch::dz
@ dz
Definition: ERF_AdvanceWDM6.cpp:272

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

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

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
1968 {
1969  // Peel back the level
1970  auto& averages = m_averages[lev];
1971 
1972  // Extrapolation is performed directly in each average's own index space;
1973  // the source fields use different centering and are not needed here.
1974 
1975  IntVect ng = averages[iavg]->nGrowVect();
1976  // Only the selected wall-normal direction is collapsed. For x/y faces
1977  // the z ghosts are tangential and must be extrapolated/fillable.
1978  ng[m_face.coordDir()] = 0;
1979 
1980  // Everything we hold was computed above
1981  if (ng.allLE(ng_fill)) { return; }
1982 
1983  IntVect ng_computed = ng_fill;
1984  ng_computed[m_face.coordDir()] = 0;
1985 
1986 #ifdef _OPENMP
1987 #pragma omp parallel if (Gpu::notInLaunchRegion())
1988 #endif
1989  for (MFIter mfi(*averages[iavg], false); mfi.isValid(); ++mfi) {
1990  // Use the average's index space. The source field and the average
1991  // have different layouts for U/V/W and for the selected wall face.
1992  const Box avg_fab = averages[iavg]->fabbox(mfi.index());
1993  Box gpbx = averages[iavg]->boxArray()[mfi.index()];
1994  gpbx.grow(ng);
1995  gpbx &= avg_fab;
1996 
1997  // Region of this average FAB that was computed from field data.
1998  Box cbx = averages[iavg]->boxArray()[mfi.index()];
1999  cbx.grow(ng_computed);
2000  cbx &= avg_fab;
2001 
2002  if (cbx.contains(gpbx)) { continue; }
2003 
2004  auto ma_arr = averages[iavg]->array(mfi);
2005 
2006  int i_lo = cbx.smallEnd(0); int i_hi = cbx.bigEnd(0);
2007  int j_lo = cbx.smallEnd(1); int j_hi = cbx.bigEnd(1);
2008  int k_lo = cbx.smallEnd(2); int k_hi = cbx.bigEnd(2);
2009  ParallelFor(gpbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
2010  {
2011  int li, lj, lk;
2012  li = i < i_lo ? i_lo : i;
2013  li = li > i_hi ? i_hi : li;
2014  lj = j < j_lo ? j_lo : j;
2015  lj = lj > j_hi ? j_hi : lj;
2016  lk = k < k_lo ? k_lo : k;
2017  lk = lk > k_hi ? k_hi : lk;
2018 
2019  ma_arr(i,j,k) = ma_arr(li,lj,lk);
2020  });
2021  } // MFiter
2022 }

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
2039 {
2040  const int dir = m_face.coordDir();
2041  const Periodicity period = tangential_periodicity(m_geom[lev], dir);
2042 
2043  if (m_terrain_type == TerrainType::EB) {
2044  mf.FillBoundary(mf.nGrowVect(), period);
2045  } else {
2046  m_planar_bndry[lev].fill(mf, period);
2047  }
2048 }
MOSTAverage::m_planar_bndry
amrex::Vector< PlanarBoundary > m_planar_bndry
Definition: ERF_MOSTAverage.H:567

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

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:559

Referenced by SurfaceLayer::get_zref().

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

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

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
733 {
734  ParmParse pp(m_pp_prefix);
735  // See zref_sentinel in ERF_MOSTAverage.H: probe with the sentinel and test the
736  // value, since queryAdd's return value stops meaning "user specified" after the
737  // first parse of the key and this routine runs once per level.
738  Real zref_tmp = zref_sentinel;
739  pp.queryAdd("most.zref",zref_tmp);
740  bool read_z = (zref_tmp > Real(0));
741  if (!read_z) { zref_tmp = zref_default; }
742  auto read_k = pp.queryarr("most.k_arr_in",m_k_in);
743 
744  const int dir = m_face.coordDir();
745  const bool is_lo_face = m_face.isLow();
746  const bool zlo = (dir == 2 && is_lo_face);
747 
748  const bool stretched = (dir == 2 && m_mesh_type != MeshType::ConstantDz);
749  if (stretched) {
750  for (int ilev : {0, lev}) {
751  const int nz = m_geom[ilev].Domain().length(2);
752  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
753  ilev < static_cast<int>(m_zlevels_stag.size()) &&
754  static_cast<int>(m_zlevels_stag[ilev].size()) == nz + 1,
755  "MOSTAverage: z levels are needed for a stretched mesh without terrain!");
756  }
757  }
758 
759  // Default behavior is to use the first cell center
760  if (!read_z && !read_k) {
761  if (stretched) {
762  const auto& zlevels = m_zlevels_stag[0];
763  zref_tmp = zlo ? cell_center_height(zlevels, 0)
764  : zlevels.back() - cell_center_height(zlevels, zlevels.size()-2);
765  } else {
766  const Real dz = m_geom[0].CellSize(dir);
767  zref_tmp = zlo ? m_geom[0].ProbLo(dir) + myhalf * dz : myhalf * dz;
768  }
769  Print() << "Reference height for MOST set to " << zref_tmp << std::endl;
770  read_z = true;
771  }
772 
773  // Specify z_ref & compute k_indx (z_ref takes precedence)
774  if (read_z && stretched) {
775  const auto& zlevels = m_zlevels_stag[lev];
776  const int nz = static_cast<int>(zlevels.size()) - 1;
777  const Real target = zlo ? zref_tmp : zlevels.back() - zref_tmp;
778  const Real first_center = cell_center_height(zlevels, 0);
779  const Real last_center = cell_center_height(zlevels, nz - 1);
780  const Real target_min = zlo ? first_center : zlevels.back() - last_center;
781  const Real target_max = zlo ? last_center : zlevels.back() - first_center;
782  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(zref_tmp >= target_min && zref_tmp <= target_max,
783  "Query point must remain inside the stretched mesh!");
784  const int lk = k_index_below(zlevels, target);
785  AMREX_ALWAYS_ASSERT_WITH_MESSAGE((zlo ? lk : nz - 1 - lk) >= m_radius,
786  "K index must be larger than averaging radius!");
787  m_k_indx[lev]->setVal(lk);
788  m_zref[lev]->setVal(zlo ? cell_center_height(zlevels, lk)
789  : zlevels.back() - cell_center_height(zlevels, lk));
790  } else if (read_z) {
791  const Real m_zlo = m_geom[lev].ProbLo(dir);
792  const Real m_zhi = m_geom[lev].ProbHi(dir);
793  const Real m_dz = m_geom[lev].CellSize(dir);
794 
795  const int dom_lo = m_geom[lev].Domain().smallEnd(dir);
796  const int dom_hi = m_geom[lev].Domain().bigEnd(dir);
797 
798  if (zlo) {
799  // Existing behavior: zlo input is an absolute coordinate.
800  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(zref_tmp >= m_zlo + myhalf * m_dz,
801  "Query point must be past first z-cell!");
802  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(zref_tmp <= m_zhi - myhalf * m_dz,
803  "Query point must be below the last z-cell!");
804 
805  const int wall_offset =
806  static_cast<int>(floor((zref_tmp - m_zlo) / m_dz - myhalf));
807 
808  AMREX_ALWAYS_ASSERT(wall_offset >= m_radius);
809 
810  m_k_indx[lev]->setVal(dom_lo + wall_offset);
811  m_zref[lev]->setVal(m_zlo +
812  (static_cast<Real>(wall_offset) + myhalf) * m_dz);
813  } else {
814  const Real domain_length = m_zhi - m_zlo;
815 
816  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(zref_tmp >= myhalf * m_dz,
817  "Query distance must reach the first cell center!");
818  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(zref_tmp <= domain_length - myhalf * m_dz,
819  "Query distance must remain inside the domain!");
820 
821  const int wall_offset = static_cast<int>(floor(zref_tmp / m_dz - myhalf));
822  const int lk = is_lo_face ? dom_lo + wall_offset : dom_hi - wall_offset;
823 
824  AMREX_ALWAYS_ASSERT(wall_offset >= m_radius);
825 
826  const Real zref_abs = is_lo_face
827  ? m_zlo + (static_cast<Real>(wall_offset) + myhalf) * m_dz
828  : m_zhi - (static_cast<Real>(wall_offset) + myhalf) * m_dz;
829  const Real zref_height = is_lo_face
830  ? zref_abs - m_zlo
831  : m_zhi - zref_abs;
832 
833  m_k_indx[lev]->setVal(lk);
834  m_zref[lev]->setVal(zref_height);
835  }
836  // Specified k_indx & compute z_ref
837  } else if (read_k) {
838  const int dom_lo = m_geom[lev].Domain().smallEnd(dir);
839  const int dom_hi = m_geom[lev].Domain().bigEnd(dir);
840  const int ncell = m_geom[lev].Domain().length(dir);
841 
842  // Preserve the established z-low convention: k_arr_in is an
843  // absolute cell index on zlo. All other faces interpret k_arr_in as
844  // a distance in cells from the selected wall.
845  const int wall_offset = zlo
846  ? m_k_in[lev] - dom_lo
847  : m_k_in[lev];
848  if (zlo) {
849  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
850  m_k_in[lev] >= dom_lo && m_k_in[lev] <= dom_hi,
851  "MOST zlo reference index must lie inside the domain!");
852  } else {
853  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
854  wall_offset >= 0 && wall_offset < ncell,
855  "MOST wall-relative reference offset must lie inside the domain!");
856  }
857  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(wall_offset >= m_radius,
858  "K index must be larger than averaging radius!");
859  const int ref_index = zlo
860  ? m_k_in[lev]
861  : (is_lo_face ? dom_lo + wall_offset : dom_hi - wall_offset);
862  m_k_indx[lev]->setVal(ref_index);
863 
864  // TODO: check that z_ref is constant across levels
865  if (dir == 2) {
866  const int dom_lo0 = m_geom[0].Domain().smallEnd(2);
867  const int dom_hi0 = m_geom[0].Domain().bigEnd(2);
868  const int wall_offset0 = zlo ? m_k_in[0] - dom_lo0 : m_k_in[0];
869  const int ref_index0 = zlo ? m_k_in[0] : dom_hi0 - wall_offset0;
870  const int k0 = ref_index0 - dom_lo0;
871  AMREX_ALWAYS_ASSERT(wall_offset0 >= m_radius);
872 
873  if (stretched) {
874  const auto& zlevels = m_zlevels_stag[0];
875  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
876  k0 >= 0 && k0 < m_geom[0].Domain().length(2),
877  "MOST reference index must lie inside the stretched domain!");
878  const Real zcell = cell_center_height(zlevels, k0);
879  m_zref[lev]->setVal(zlo ? zcell : zlevels.back() - zcell);
880  } else {
881  const Real dz0 = m_geom[0].CellSize(2);
882  const Real zlo0 = m_geom[0].ProbLo(2);
883  const Real zhi0 = m_geom[0].ProbHi(2);
884  const Real zref_abs = is_lo_face
885  ? zlo0 + (static_cast<Real>(wall_offset0) + myhalf) * dz0
886  : zhi0 - (static_cast<Real>(wall_offset0) + myhalf) * dz0;
887  m_zref[lev]->setVal(zlo
888  ? zref_abs
889  : (is_lo_face ? zref_abs - zlo0 : zhi0 - zref_abs));
890  }
891  } else {
892  const Real m_dz = m_geom[lev].CellSize(dir);
893  const Real m_zlo = m_geom[lev].ProbLo(dir);
894  const Real m_zhi = m_geom[lev].ProbHi(dir);
895  const Real zref_abs = is_lo_face
896  ? m_zlo + (static_cast<Real>(wall_offset) + myhalf) * m_dz
897  : m_zhi - (static_cast<Real>(wall_offset) + myhalf) * m_dz;
898 
899  m_zref[lev]->setVal(zlo
900  ? zref_abs
901  : (is_lo_face ? zref_abs - m_zlo : m_zhi - zref_abs));
902  }
903  }
904 }
length
const Real length
Definition: ERF_InitCustomPert_AnelasticWallDiffusion.H:14
pp
ParmParse pp("prob")
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:615
MOSTAverage::m_pp_prefix
std::string m_pp_prefix
Definition: ERF_MOSTAverage.H:547
MOSTAverage::m_mesh_type
MeshType m_mesh_type
Definition: ERF_MOSTAverage.H:548
MOSTAverage::m_zlevels_stag
amrex::Vector< amrex::Vector< amrex::Real > > m_zlevels_stag
Definition: ERF_MOSTAverage.H:550
MOSTAverage::m_k_in
amrex::Vector< int > m_k_in
Definition: ERF_MOSTAverage.H:579

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

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
1142 {
1143  // Peel back the level
1144  auto& fields = m_fields[lev];
1145 
1146  // MFIter over CC data
1147  int imf_cc = 3;
1148 
1149  ParmParse pp(m_pp_prefix);
1150  // See zref_sentinel in ERF_MOSTAverage.H: probe with the sentinel and test the
1151  // value, since queryAdd's return value stops meaning "user specified" after the
1152  // first parse of the key and this routine runs once per level.
1153  Real zref_tmp = zref_sentinel;
1154  pp.queryAdd("most.zref",zref_tmp);
1155  bool read_zref = (zref_tmp > Real(0));
1156  if (!read_zref) {
1157  zref_tmp = zref_default;
1158  Print() << "most.zref not specified, query distance default is " << zref_tmp << std::endl;
1159  }
1160 
1161  AMREX_ALWAYS_ASSERT(m_face.coordDir() == 2);
1162  int klo = m_geom[lev].Domain().smallEnd(2);
1163 
1164  // Capture for device
1165  Real d_zref = zref_tmp;
1166  Real d_radius = static_cast<Real>(m_radius);
1167 
1168  const auto dxInv = m_geom[lev].InvCellSizeArray();
1169 
1170  // The indices are needed everywhere an average is computed, ghost cells
1171  // included; the box is made nodal so that we also cover the U & V averages,
1172  // which are face centered
1173  const IntVect ng_indx = max(get_ng_fill(lev), IntVect(1,1,0));
1174 
1175  for (MFIter mfi(*fields[imf_cc], TileNoZ()); mfi.isValid(); ++mfi) {
1176  Box npbx = mfi.tilebox(IntVect(1,1,0),IntVect(1,1,0));
1177  const Box vbx = mfi.validbox();
1178  npbx = mfi.tilebox(IntVect(1,1,0),ng_indx);
1179 
1180  if (vbx.smallEnd(2) != klo) { continue; }
1181 
1182  int kmax = npbx.bigEnd(2);
1183 
1184  npbx.makeSlab(2,klo);
1185 
1186  const auto z_phys_arr = m_z_phys_nd[lev]->const_array(mfi);
1187  auto i_arr = m_i_indx[lev]->array(mfi);
1188  auto j_arr = m_j_indx[lev]->array(mfi);
1189  auto k_arr = m_k_indx[lev]->array(mfi);
1190  auto zref_arr = m_zref[lev]->array(mfi);
1191  ParallelFor(npbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
1192  {
1193  // Elements of normal vector
1194  Real met_h_xi = Compute_h_xi_AtCellCenter (i,j,k,dxInv,z_phys_arr);
1195  Real met_h_eta = Compute_h_eta_AtCellCenter(i,j,k,dxInv,z_phys_arr);
1196  Real mag = std::sqrt(met_h_xi*met_h_xi + met_h_eta*met_h_eta + one);
1197 
1198  // Unit-normal vector scaled by z_ref
1199  Real delta_x = -met_h_xi/mag * d_zref;
1200  Real delta_y = -met_h_eta/mag * d_zref;
1201  Real delta_z = one/mag * d_zref;
1202 
1203  // Compute i & j as displacements (no grid stretching)
1204  int delta_i = static_cast<int>(std::round(delta_x*dxInv[0]));
1205  int delta_j = static_cast<int>(std::round(delta_y*dxInv[1]));
1206  int i_new = i + delta_i;
1207  int j_new = j + delta_j;
1208  i_arr(i,j,0) = i_new;
1209  j_arr(i,j,0) = j_new;
1210 
1211  // Search for k (grid is stretched in z)
1212  Real z_bot_face = fourth * ( z_phys_arr(i ,j ,k) + z_phys_arr(i+1,j ,k)
1213  + z_phys_arr(i ,j+1,k) + z_phys_arr(i+1,j+1,k) );
1214  Real z_target = z_bot_face + delta_z;
1215  k_arr(i,j,0) = klo;
1216  zref_arr(i,j,0) = myhalf * z_bot_face +
1217  Real(0.125) * ( z_phys_arr(i ,j ,k+1) + z_phys_arr(i+1,j ,k+1)
1218  + z_phys_arr(i ,j+1,k+1) + z_phys_arr(i+1,j+1,k+1) );
1219  for (int lk(klo); lk<=kmax; ++lk) {
1220  Real z_lo = fourth * ( z_phys_arr(i_new,j_new ,lk ) + z_phys_arr(i_new+1,j_new ,lk )
1221  + z_phys_arr(i_new,j_new+1,lk ) + z_phys_arr(i_new+1,j_new+1,lk ) );
1222  Real z_hi = fourth * ( z_phys_arr(i_new,j_new ,lk+1) + z_phys_arr(i_new+1,j_new ,lk+1)
1223  + z_phys_arr(i_new,j_new+1,lk+1) + z_phys_arr(i_new+1,j_new+1,lk+1) );
1224  if (in_cell_z(z_target, z_lo, z_hi)) {
1225  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(lk >= d_radius,
1226  "K index must be larger than averaging radius!");
1227  amrex::ignore_unused(d_radius);
1228  k_arr(i,j,0) = lk;
1229  zref_arr(i,j,0) = myhalf * (z_hi + z_lo) - z_bot_face;
1230  break;
1231  }
1232  }
1233  });
1234  }
1235 }
klo
const int klo
Definition: ERF_InitCustomPert_ABL.H:75
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:215
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:240
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
1325 {
1326  AMREX_ALWAYS_ASSERT(m_face.coordDir() == 2);
1327 
1328  // Peel back the level
1329  auto& fields = m_fields[lev];
1330 
1331  // MFIter over CC data
1332  int imf_cc = 3;
1333 
1334  ParmParse pp(m_pp_prefix);
1335  // See zref_sentinel in ERF_MOSTAverage.H: probe with the sentinel and test the
1336  // value, since queryAdd's return value stops meaning "user specified" after the
1337  // first parse of the key and this routine runs once per level.
1338  Real zref_tmp = zref_sentinel;
1339  pp.queryAdd("most.zref",zref_tmp);
1340  bool read_zref = (zref_tmp > Real(0));
1341  if (!read_zref) {
1342  zref_tmp = zref_default;
1343  Print() << "most.zref not specified, query distance default is " << zref_tmp << std::endl;
1344  }
1345  int klo = m_geom[lev].Domain().smallEnd(2);
1346 
1347  // Capture for device
1348  Real d_zref = zref_tmp;
1349  const auto plo = m_geom[lev].ProbLoArray();
1350 
1351  RealVect base;
1352  const auto dx = m_geom[lev].CellSizeArray();
1353  const auto dxInv = m_geom[lev].InvCellSizeArray();
1354  IntVect ng = m_x_pos[lev]->nGrowVect(); ng[2]=0;
1355  const int position_ng = (m_radius > 1) ? m_radius : 1;
1356 
1357  // The positions are read over the averaging stencil at every cell where an
1358  // average is computed, ghost cells included; the box is made nodal so that
1359  // we also cover the U & V averages, which are face centered
1360  const IntVect ng_pos = max(get_ng_fill(lev) + IntVect(m_radius,m_radius,0),
1361  IntVect(position_ng,position_ng,0));
1362 
1363  for (MFIter mfi(*fields[imf_cc], TileNoZ()); mfi.isValid(); ++mfi) {
1364  Box npbx = mfi.tilebox(IntVect(1,1,0),ng_pos);
1365  Box gtbx = mfi.growntilebox(ng);
1366  RealBox grb{gtbx,dx.data(),base.dataPtr()};
1367 
1368  if (mfi.validbox().smallEnd(2) != klo) { continue; }
1369 
1370  npbx.makeSlab(2,klo);
1371 
1372  const auto z_phys_arr = m_z_phys_nd[lev]->const_array(mfi);
1373  auto x_pos_arr = m_x_pos[lev]->array(mfi);
1374  auto y_pos_arr = m_y_pos[lev]->array(mfi);
1375  auto z_pos_arr = m_z_pos[lev]->array(mfi);
1376  auto zref_arr = m_zref[lev]->array(mfi);
1377  ParallelFor(npbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
1378  {
1379  // Elements of normal vector
1380  Real met_h_xi = Compute_h_xi_AtCellCenter (i,j,k,dxInv,z_phys_arr);
1381  Real met_h_eta = Compute_h_eta_AtCellCenter(i,j,k,dxInv,z_phys_arr);
1382  Real imag = one / std::sqrt(met_h_xi*met_h_xi + met_h_eta*met_h_eta + one);
1383 
1384  // Unit-normal vector scaled by z_ref
1385  Real delta_x = -met_h_xi * imag * d_zref;
1386  Real delta_y = -met_h_eta * imag * d_zref;
1387  Real delta_z = imag * d_zref;
1388 
1389  // Position of the current node (indx:0,0,1)
1390  Real x0 = plo[0] + ((Real) i + myhalf) * dx[0];
1391  Real y0 = plo[1] + ((Real) j + myhalf) * dx[1];
1392 
1393  // Final position at end of vector
1394  x_pos_arr(i,j,0) = x0 + delta_x;
1395  y_pos_arr(i,j,0) = y0 + delta_y;
1396  Real z_bot_face = fourth * ( z_phys_arr(i ,j ,k) + z_phys_arr(i+1,j ,k)
1397  + z_phys_arr(i ,j+1,k) + z_phys_arr(i+1,j+1,k) );
1398  z_pos_arr(i,j,0) = z_bot_face + delta_z;
1399 
1400  // NOTE: Normal vector end point can be below the surface for concave regions.
1401  // Here we protect against that by augmenting the normal if needed.
1402  int i_new = (int) ((x_pos_arr(i,j,0) - plo[0]) / dx[0] - myhalf);
1403  int j_new = (int) ((y_pos_arr(i,j,0) - plo[1]) / dx[1] - myhalf);
1404  Real z_new_bot_face = fourth * ( z_phys_arr(i_new,j_new ,k) + z_phys_arr(i_new+1,j_new ,k)
1405  + z_phys_arr(i_new,j_new+1,k) + z_phys_arr(i_new+1,j_new+1,k) );
1406  if (z_pos_arr(i,j,0) < z_new_bot_face) {
1407  z_pos_arr(i,j,0) = z_new_bot_face + delta_z;
1408  }
1409 
1410  zref_arr(i,j,0) = delta_z;
1411 
1412  // Destination position must be contained on the current process!
1413  Real pos[] = {x_pos_arr(i,j,0)-plo[0],y_pos_arr(i,j,0)-plo[1],myhalf*dx[2]};
1414  amrex::ignore_unused(pos);
1415  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(grb.contains(&pos[0]),
1416  "Query point outside of proc domain!");
1417  });
1418  }
1419 }

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

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

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
959 {
960  // See zref_sentinel in ERF_MOSTAverage.H: probe with the sentinel and test the
961  // value, since queryAdd's return value stops meaning "user specified" after the
962  // first parse of the key and this routine runs once per level.
963  Real zref_tmp = zref_sentinel;
964  ParmParse pp(m_pp_prefix);
965  pp.queryAdd("most.zref",zref_tmp);
966  bool read_z = (zref_tmp > Real(0));
967  if (!read_z) { zref_tmp = zref_default; }
968 
969  if (read_z) {
970  m_zref[lev]->setVal( zref_tmp );
971  // Default behavior is to use the first cell center
972  } else {
973  Real m_dz = m_geom[0].CellSize(2);
974  zref_tmp = myhalf * m_dz;
975  m_zref[lev]->setVal( zref_tmp );
976  Print() << "Reference height for MOST set to " << zref_tmp << std::endl;
977  }
978 }

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
1245 {
1246  AMREX_ALWAYS_ASSERT(m_face.coordDir() == 2);
1247 
1248  // Peel back the level
1249  auto& fields = m_fields[lev];
1250 
1251  // MFIter over CC data
1252  int imf_cc = 3;
1253 
1254  ParmParse pp(m_pp_prefix);
1255  // See zref_sentinel in ERF_MOSTAverage.H: probe with the sentinel and test the
1256  // value, since queryAdd's return value stops meaning "user specified" after the
1257  // first parse of the key and this routine runs once per level.
1258  Real zref_tmp = zref_sentinel;
1259  pp.queryAdd("most.zref",zref_tmp);
1260  bool read_zref = (zref_tmp > Real(0));
1261  if (!read_zref) {
1262  zref_tmp = zref_default;
1263  Print() << "most.zref not specified, query distance default is " << zref_tmp << std::endl;
1264  } else {
1265  m_zref[lev]->setVal(zref_tmp);
1266  }
1267  int klo = m_geom[lev].Domain().smallEnd(2);
1268 
1269  // Capture for device
1270  Real d_zref = zref_tmp;
1271  const auto plo = m_geom[lev].ProbLoArray();
1272 
1273  RealVect base;
1274  const auto dx = m_geom[lev].CellSizeArray();
1275  IntVect ng = m_x_pos[lev]->nGrowVect(); ng[2]=0;
1276  const int position_ng = (m_radius > 1) ? m_radius : 1;
1277 
1278  // The positions are read over the averaging stencil at every cell where an
1279  // average is computed, ghost cells included; the box is made nodal so that
1280  // we also cover the U & V averages, which are face centered
1281  const IntVect ng_pos = max(get_ng_fill(lev) + IntVect(m_radius,m_radius,0),
1282  IntVect(position_ng,position_ng,0));
1283 
1284  for (MFIter mfi(*fields[imf_cc], TileNoZ()); mfi.isValid(); ++mfi) {
1285  Box npbx = mfi.tilebox(IntVect(1,1,0),ng_pos);
1286  const Box vbx = mfi.validbox();
1287  Box gtbx = mfi.growntilebox(ng);
1288 
1289  if (vbx.smallEnd(2) != klo) { continue; }
1290 
1291  npbx.makeSlab(2,klo);
1292 
1293  RealBox grb{gtbx,dx.data(),base.dataPtr()};
1294 
1295  const auto z_phys_arr = m_z_phys_nd[lev]->const_array(mfi);
1296  auto x_pos_arr = m_x_pos[lev]->array(mfi);
1297  auto y_pos_arr = m_y_pos[lev]->array(mfi);
1298  auto z_pos_arr = m_z_pos[lev]->array(mfi);
1299  ParallelFor(npbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
1300  {
1301  // Final position at end of vector
1302  x_pos_arr(i,j,k) = plo[0] + ((Real) i + myhalf) * dx[0];
1303  y_pos_arr(i,j,k) = plo[1] + ((Real) j + myhalf) * dx[1];
1304  Real z_bot_face = fourth * ( z_phys_arr(i ,j ,k) + z_phys_arr(i+1,j ,k)
1305  + z_phys_arr(i ,j+1,k) + z_phys_arr(i+1,j+1,k) );
1306  z_pos_arr(i,j,k) = z_bot_face + d_zref;
1307 
1308  // Destination position must be contained on the current process!
1309  Real pos[] = {x_pos_arr(i,j,k)-plo[0],y_pos_arr(i,j,k)-plo[1],myhalf*dx[2]};
1310  amrex::ignore_unused(pos);
1311  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(grb.contains(&pos[0]),
1312  "Query point outside of proc domain!");
1313  });
1314  }
1315 }

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

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
446 {
447  m_fields[lev][0] = &vars_old[lev][Vars::xvel];
448  m_fields[lev][1] = &vars_old[lev][Vars::yvel];
449  m_fields[lev][2] = &vars_old[lev][Vars::zvel];
450  m_fields[lev][3] = Theta_prim[lev].get();
451  m_fields[lev][4] = Qv_prim[lev].get();
452  m_fields[lev][5] = Qr_prim[lev].get();
453 }

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.

920 {
921  const int dir = m_face.coordDir();
922  if (dir == 2 || m_interp || m_norm_vec || !m_k_indx[lev]) { return; }
923 
924  const MultiFab& source = *m_fields[lev][3];
925  const int face_index = m_face.isLow()
926  ? m_geom[lev].Domain().smallEnd(dir)
927  : m_geom[lev].Domain().bigEnd(dir);
928  const int reference = m_k_indx[lev]->min(0);
929 
930  for (MFIter mfi(source, false); mfi.isValid(); ++mfi) {
931  const Box& valid_box = mfi.validbox();
932  const bool owns_face = m_face.isLow()
933  ? valid_box.smallEnd(dir) == face_index
934  : valid_box.bigEnd(dir) == face_index;
935  if (!owns_face) { continue; }
936 
937  Box required = valid_box;
938  required.setSmall(dir, reference - m_radius);
939  required.setBig(dir, reference + m_radius);
940  required.growHi(0, 1);
941  required.growHi(1, 1);
942  required.growHi(2, 1);
943 
944  Box available = valid_box;
945  available.grow(source.nGrowVect());
946  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
947  available.contains(required),
948  "MOST lateral reference plus averaging radius is outside the face-owned FAB.");
949  }
950 }

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
3191 {
3192  // Peel back the level
3193  auto& fields = m_fields[lev];
3194  auto& averages = m_averages[lev];
3195 
3196  // MFIter on CC
3197  int imf_cc = 3;
3198 
3199 
3200  int navg = m_navg - 1;
3201 
3202  std::ofstream ofile;
3203  ofile.open (std::string("MOST_averages_" + std::to_string(m_face) + ".txt"));
3204  ofile << "Averages computed via MOSTAverages class:\n";
3205 
3206  const int dir = m_face.coordDir();
3207  int sm_index;
3208  if (m_face.isLow()) {
3209  sm_index = m_geom[lev].Domain().smallEnd(dir);
3210  } else {
3211  sm_index = m_geom[lev].Domain().bigEnd(dir);
3212  }
3213 
3214  for (MFIter mfi(*fields[imf_cc], TileNoZ()); mfi.isValid(); ++mfi) {
3215  Box pbx = mfi.tilebox();
3216 
3217  if (m_face.isLow()) {
3218  if (mfi.validbox().smallEnd(dir) != sm_index ||
3219  pbx.smallEnd(dir) != sm_index) {
3220  continue;
3221  }
3222  } else {
3223  if (mfi.validbox().bigEnd(dir) != sm_index ||
3224  pbx.bigEnd(dir) != sm_index) {
3225  continue;
3226  }
3227  }
3228 
3229  pbx.setSmall(dir, sm_index); pbx.setBig(dir, sm_index);
3230 
3231  int il = pbx.smallEnd(0); int iu = pbx.bigEnd(0);
3232  int jl = pbx.smallEnd(1); int ju = pbx.bigEnd(1);
3233  int kl = pbx.smallEnd(2); int ku = pbx.bigEnd(2);
3234 
3235  for (int k(kl); k <= ku; ++k) {
3236  for (int j(jl); j <= ju; ++j) {
3237  for (int i(il); i <= iu; ++i) {
3238  if (dir == 0) {
3239  ofile << "(J,K): " << "(" << j << "," << k << ")" << "\n";
3240  } else if (dir == 1) {
3241  ofile << "(I,K): " << "(" << i << "," << k << ")" << "\n";
3242  } else {
3243  ofile << "(I,J): " << "(" << i << "," << j << ")" << "\n";
3244  }
3245 
3246  for (int iavg(0); iavg <= navg; ++iavg) {
3247  auto mf_arr = averages[iavg]->array(mfi);
3248  ofile << "iavg val: "
3249  << iavg << ' '
3250  << mf_arr(i,j,k) << "\n";
3251  }
3252  ofile << "\n";
3253  }
3254  }
3255  }
3256  }
3257  ofile.close();
3258 }
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
2993 {
2994  // Peel back the level
2995  auto& fields = m_fields[lev];
2996  auto& k_indx = m_k_indx[lev];
2997 
2998  // MFIter on CC
2999  int imf_cc = 3;
3000 
3001  std::ofstream ofile;
3002  ofile.open ("MOST_k_indices_" + std::to_string(m_face) + ".txt");
3003  ofile << "K indices used to compute averages via MOSTAverages class:\n";
3004 
3005  const int dir = m_face.coordDir();
3006  int sm_index;
3007  if (m_face.isLow()) {
3008  sm_index = m_geom[lev].Domain().smallEnd(dir);
3009  } else {
3010  sm_index = m_geom[lev].Domain().bigEnd(dir);
3011  }
3012 
3013  for (MFIter mfi(*fields[imf_cc], TileNoZ()); mfi.isValid(); ++mfi) {
3014  Box pbx = mfi.tilebox();
3015 
3016  if (m_face.isLow()) {
3017  if (mfi.validbox().smallEnd(dir) != sm_index ||
3018  pbx.smallEnd(dir) != sm_index) {
3019  continue;
3020  }
3021  } else {
3022  if (mfi.validbox().bigEnd(dir) != sm_index ||
3023  pbx.bigEnd(dir) != sm_index) {
3024  continue;
3025  }
3026  }
3027 
3028  pbx.setSmall(dir, sm_index); pbx.setBig(dir, sm_index);
3029 
3030  int il = pbx.smallEnd(0); int iu = pbx.bigEnd(0);
3031  int jl = pbx.smallEnd(1); int ju = pbx.bigEnd(1);
3032  int kl = pbx.smallEnd(2); int ku = pbx.bigEnd(2);
3033 
3034  auto k_arr = k_indx->array(mfi);
3035  for (int k(kl); k <= ku; ++k) {
3036  for (int j(jl); j <= ju; ++j) {
3037  for (int i(il); i <= iu; ++i) {
3038  if (dir == 0) {
3039  ofile << "(J,K): " << "(" << j << "," << k << ")" << "\n";
3040  } else if (dir == 1) {
3041  ofile << "(I,K): " << "(" << i << "," << k << ")" << "\n";
3042  } else {
3043  ofile << "(I,J): " << "(" << i << "," << j << ")" << "\n";
3044  }
3045  ofile << "K_ind: "
3046  << k_arr(i,j,k) << "\n";
3047  ofile << "\n";
3048  }
3049  }
3050  }
3051  }
3052  ofile.close();
3053 }
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
3063 {
3064  // Peel back the level
3065  auto& fields = m_fields[lev];
3066  auto& k_indx = m_k_indx[lev];
3067  auto& j_indx = m_j_indx[lev];
3068  auto& i_indx = m_i_indx[lev];
3069 
3070  // MFIter on CC
3071  int imf_cc = 3;
3072 
3073  std::ofstream ofile;
3074  ofile.open ("MOST_ijk_indices_" + std::to_string(m_face) + ".txt");
3075  ofile << "IJK indices used to compute averages via MOSTAverages class:\n";
3076 
3077  const int dir = m_face.coordDir();
3078  int sm_index;
3079  if (m_face.isLow()) {
3080  sm_index = m_geom[lev].Domain().smallEnd(dir);
3081  } else {
3082  sm_index = m_geom[lev].Domain().bigEnd(dir);
3083  }
3084 
3085  for (MFIter mfi(*fields[imf_cc], TileNoZ()); mfi.isValid(); ++mfi) {
3086  Box pbx = mfi.tilebox();
3087 
3088  if (m_face.isLow()) {
3089  if (mfi.validbox().smallEnd(dir) != sm_index ||
3090  pbx.smallEnd(dir) != sm_index) {
3091  continue;
3092  }
3093  } else {
3094  if (mfi.validbox().bigEnd(dir) != sm_index ||
3095  pbx.bigEnd(dir) != sm_index) {
3096  continue;
3097  }
3098  }
3099 
3100  pbx.setSmall(dir, sm_index); pbx.setBig(dir, sm_index);
3101 
3102  int il = pbx.smallEnd(0); int iu = pbx.bigEnd(0);
3103  int jl = pbx.smallEnd(1); int ju = pbx.bigEnd(1);
3104  int kl = pbx.smallEnd(2); int ku = pbx.bigEnd(2);
3105 
3106  auto k_arr = k_indx->array(mfi);
3107  auto j_arr = j_indx ? j_indx->array(mfi) : Array4<int> {};
3108  auto i_arr = i_indx ? i_indx->array(mfi) : Array4<int> {};
3109 
3110  for (int k(kl); k <= ku; ++k) {
3111  for (int j(jl); j <= ju; ++j) {
3112  for (int i(il); i <= iu; ++i) {
3113  ofile << "(I1,J1,K1): " << "(" << i << "," << j << "," << k << ")" << "\n";
3114 
3115  // This diagnostic reports the cell-centered reference used
3116  // for the CC fields. The high-face +1 adjustment applies
3117  // only to the face-centered velocity normal to the wall.
3118  const int ref = k_arr(i,j,k);
3119  int im = i_arr ? i_arr(i,j,k) : i;
3120  int jm = j_arr ? j_arr(i,j,k) : j;
3121  int km = k;
3122  if (dir == 0) {
3123  im = ref;
3124  } else if (dir == 1) {
3125  jm = ref;
3126  } else {
3127  km = ref;
3128  }
3129 
3130  ofile << "(I2,J2,K2): "
3131  << "(" << im << "," << jm << "," << km << ")" << "\n";
3132  ofile << "\n";
3133  }
3134  }
3135  }
3136  }
3137  ofile.close();
3138 }
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
3150 {
3151  // Peel back the level
3152  auto& fields = m_fields[lev];
3153  auto& x_pos_mf = m_x_pos[lev];
3154  auto& z_pos_mf = m_z_pos[lev];
3155 
3156  // MFIter on CC
3157  int imf_cc = 3;
3158 
3159  int klo = m_geom[lev].Domain().smallEnd(2);
3160 
3161  std::ofstream ofile;
3162  ofile.open ("MOST_xz_positions.txt");
3163 
3164  for (MFIter mfi(*fields[imf_cc], TileNoZ()); mfi.isValid(); ++mfi) {
3165  Box pbx = mfi.tilebox();
3166 
3167  if(pbx.smallEnd(2) != klo) { continue; }
3168 
3169  pbx.makeSlab(2,klo);
3170 
3171  int il = pbx.smallEnd(0); int iu = pbx.bigEnd(0);
3172 
3173  auto x_pos_arr = x_pos_mf->array(mfi);
3174  auto z_pos_arr = z_pos_mf->array(mfi);
3175 
3176  int k = 0;
3177  for (int i(il); i <= iu; ++i)
3178  ofile << x_pos_arr(i,j,k) << ' ' << z_pos_arr(i,j,k) << "\n";
3179  }
3180  ofile.close();
3181 }
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: