ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
Public Member Functions | Static Public Member Functions | Protected Attributes | List of all members
MOSTAverage Class Reference

#include <ERF_MOSTAverage.H>

Collaboration diagram for MOSTAverage:
Collaboration graph
[legend]

Public Member Functions

 MOSTAverage (amrex::Vector< amrex::Geometry > geom, const bool &has_zphys, std::string a_pp_prefix, const TerrainType &m_terrain_type)
 
 ~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_region_normalization (const int &)
 
void set_k_indices_N (const int &lev)
 
void set_k_indices_T (const int &lev)
 
void set_norm_indices_T (const int &lev)
 
void set_z_positions_T (const int &lev)
 
void set_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 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::Real get_zref () const
 

Static Public Member Functions

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

Protected Attributes

const amrex::Vector< amrex::Geometry > m_geom
 
amrex::Vector< amrex::Vector< amrex::MultiFab * > > m_fields
 
amrex::Vector< amrex::MultiFab * > m_z_phys_nd
 
std::string m_pp_prefix
 
TerrainType m_terrain_type
 
int m_nvar {6}
 
int m_navg {6}
 
int m_maxlev {0}
 
int m_policy {0}
 
bool m_rotate {false}
 
amrex::Real m_zref {10.0}
 
amrex::Vector< std::unique_ptr< amrex::MultiFab > > m_x_pos
 
amrex::Vector< std::unique_ptr< amrex::MultiFab > > m_y_pos
 
amrex::Vector< std::unique_ptr< amrex::MultiFab > > m_z_pos
 
amrex::Vector< std::unique_ptr< amrex::iMultiFab > > m_i_indx
 
amrex::Vector< std::unique_ptr< amrex::iMultiFab > > m_j_indx
 
amrex::Vector< std::unique_ptr< amrex::iMultiFab > > m_k_indx
 
amrex::Vector< amrex::Vector< std::unique_ptr< amrex::MultiFab > > > m_averages
 
amrex::Vector< amrex::Vector< std::unique_ptr< amrex::MultiFab > > > m_rot_fields
 
amrex::Vector< amrex::Vector< int > > m_ncell_plane
 
amrex::Vector< amrex::Vector< amrex::Real > > m_plane_average
 
int m_radius {0}
 
int m_ncell_region {1}
 
amrex::Vector< int > m_k_in
 
bool m_interp {false}
 
bool m_norm_vec {false}
 
bool m_t_avg {false}
 
amrex::Vector< int > m_t_init
 
amrex::Real 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
 

Constructor & Destructor Documentation

◆ MOSTAverage() [1/3]

MOSTAverage::MOSTAverage ( amrex::Vector< amrex::Geometry >  geom,
const bool &  has_zphys,
std::string  a_pp_prefix,
const TerrainType &  m_terrain_type 
)
explicit

◆ ~MOSTAverage()

MOSTAverage::~MOSTAverage ( )
inline
22  {}

◆ MOSTAverage() [2/3]

MOSTAverage::MOSTAverage ( MOSTAverage &&  )
defaultnoexcept

◆ MOSTAverage() [3/3]

MOSTAverage::MOSTAverage ( const MOSTAverage &  other)
delete

Member Function Documentation

◆ compute_averages()

void MOSTAverage::compute_averages ( const int &  lev)

Function to call the type of average computation.

Parameters
[in]levCurrent level
666 {
667  if (m_rotate) set_rotated_fields(lev);
668 
669  switch(m_policy) {
670  case 0: // Standard plane average
671  compute_plane_averages(lev);
672  break;
673  case 1: // Local region/point
674  compute_region_averages(lev);
675  break;
676  default:
677  AMREX_ASSERT_WITH_MESSAGE(false, "Unknown policy for MOSTAverage!");
678  }
679 
680  // We have initialized the averages
681  if (m_t_avg) m_t_init[lev] = 1;
682 }
MOSTAverage::m_t_avg
bool m_t_avg
Definition: ERF_MOSTAverage.H:218
MOSTAverage::compute_plane_averages
void compute_plane_averages(const int &lev)
Definition: ERF_MOSTAverage.cpp:691
MOSTAverage::m_policy
int m_policy
Definition: ERF_MOSTAverage.H:188
MOSTAverage::compute_region_averages
void compute_region_averages(const int &lev)
Definition: ERF_MOSTAverage.cpp:976
MOSTAverage::m_t_init
amrex::Vector< int > m_t_init
Definition: ERF_MOSTAverage.H:219
MOSTAverage::set_rotated_fields
void set_rotated_fields(const int &lev)
Definition: ERF_MOSTAverage.cpp:263
MOSTAverage::m_rotate
bool m_rotate
Definition: ERF_MOSTAverage.H:189
Here is the call graph for this function:

◆ compute_plane_averages()

void MOSTAverage::compute_plane_averages ( const int &  lev)

Function to compute average over a plane.

Parameters
[in]levCurrent level
692 {
693  // Peel back the level
694  auto& fields = m_fields[lev];
695  auto& rot_fields = m_rot_fields[lev];
696  auto& averages = m_averages[lev];
697  const auto & geom = m_geom[lev];
698 
699  auto& z_phys = m_z_phys_nd[lev];
700  auto& x_pos = m_x_pos[lev];
701  auto& y_pos = m_y_pos[lev];
702  auto& z_pos = m_z_pos[lev];
703 
704  auto& i_indx = m_i_indx[lev];
705  auto& j_indx = m_j_indx[lev];
706  auto& k_indx = m_k_indx[lev];
707 
708  auto& ncell_plane = m_ncell_plane[lev];
709  auto& plane_average = m_plane_average[lev];
710 
711  // Set factors for time averaging
712  Real d_fact_new, d_fact_old;
713  if (m_t_avg && m_t_init[lev]) {
714  d_fact_new = m_fact_new;
715  d_fact_old = m_fact_old;
716  } else {
717  d_fact_new = 1.0;
718  d_fact_old = 0.0;
719  }
720 
721  // GPU array to accumulate averages into
722  Gpu::DeviceVector<Real> pavg(plane_average.size(), 0.0);
723  Real* plane_avg = pavg.data();
724 
725  // Vectors for normalization and buffer storage
726  Vector<Real> denom(plane_average.size(),0.0);
727  Vector<Real> val_old(plane_average.size(),0.0);
728 
729  //
730  //----------------------------------------------------------
731  // Averages over all the fields
732  //----------------------------------------------------------
733  //
734  Box domain = geom.Domain();
735 
736  Array<int,AMREX_SPACEDIM> is_per = {0,0,0};
737  for (int idim(0); idim < AMREX_SPACEDIM-1; ++idim) {
738  if (geom.isPeriodic(idim)) is_per[idim] = 1;
739  }
740 
741  // Averages for U,V,T,Qv (not Qr or W)
742  for (int imf(0); imf < 4; ++imf) {
743 
744  // Continue if no valid Qv pointer
745  if (!fields[imf]) continue;
746 
747  denom[imf] = 1.0 / (Real)ncell_plane[imf];
748  val_old[imf] = plane_average[imf]*d_fact_old;
749 
750 #ifdef _OPENMP
751 #pragma omp parallel if (Gpu::notInLaunchRegion())
752 #endif
753  for (MFIter mfi(*fields[imf], TileNoZ()); mfi.isValid(); ++mfi) {
754  Box vbx = mfi.validbox(); // This is the grid (not tile)
755  Box pbx = mfi.tilebox(); // This is the tile (not grid)
756  pbx.setSmall(2,0); pbx.setBig(2,0);
757 
758  // Avoid double counting nodal data by changing the high end when we are
759  // at the high side of the grid (not just of the tile)
760  IndexType ixt = averages[imf]->boxArray().ixType();
761  for (int idim(0); idim < AMREX_SPACEDIM-1; ++idim) {
762  if ( ixt.nodeCentered(idim) && (pbx.bigEnd(idim) == vbx.bigEnd(idim)) ) {
763  int dom_hi = domain.bigEnd(idim)+1;
764  if (pbx.bigEnd(idim) < dom_hi || is_per[idim]) {
765  pbx.growHi(idim,-1);
766  }
767  }
768  }
769 
770  auto mf_arr = (m_rotate) ? rot_fields[imf]->const_array(mfi) :
771  fields[imf]->const_array(mfi);
772 
773  if (m_interp) {
774  const auto plo = geom.ProbLoArray();
775  const auto dxInv = geom.InvCellSizeArray();
776  const auto z_phys_arr = z_phys->const_array(mfi);
777  auto x_pos_arr = x_pos->array(mfi);
778  auto y_pos_arr = y_pos->array(mfi);
779  auto z_pos_arr = z_pos->array(mfi);
780  ParallelFor(Gpu::KernelInfo().setReduction(true), pbx, [=]
781  AMREX_GPU_DEVICE(int i, int j, int k, Gpu::Handler const& handler) noexcept
782  {
783  Real interp{0};
784  trilinear_interp_T(x_pos_arr(i,j,k), y_pos_arr(i,j,k), z_pos_arr(i,j,k),
785  &interp, mf_arr, z_phys_arr, plo, dxInv, 1);
786  Real val = interp;
787  Gpu::deviceReduceSum(&plane_avg[imf], val, handler);
788  });
789  } else {
790  auto k_arr = k_indx->const_array(mfi);
791  auto j_arr = j_indx ? j_indx->const_array(mfi) : Array4<const int> {};
792  auto i_arr = i_indx ? i_indx->const_array(mfi) : Array4<const int> {};
793  ParallelFor(Gpu::KernelInfo().setReduction(true), pbx, [=]
794  AMREX_GPU_DEVICE(int i, int j, int k, Gpu::Handler const& handler) noexcept
795  {
796  int mk = k_arr(i,j,k);
797  int mj = j_arr ? j_arr(i,j,k) : j;
798  int mi = i_arr ? i_arr(i,j,k) : i;
799  Real val = mf_arr(mi,mj,mk);
800  Gpu::deviceReduceSum(&plane_avg[imf], val, handler);
801  });
802  }
803  }
804  }
805 
806  //
807  //------------------------------------------------------------------------
808  // Averages for virtual potential temperature
809  // (This is cell-centered so we don't need to worry about double-counting)
810  //------------------------------------------------------------------------
811  //
812  if (fields[3]) // We have water vapor
813  {
814  int iavg = 4;
815  denom[iavg] = 1.0 / (Real)ncell_plane[iavg];
816  val_old[iavg] = plane_average[iavg]*d_fact_old;
817 
818 #ifdef _OPENMP
819 #pragma omp parallel if (Gpu::notInLaunchRegion())
820 #endif
821  for (MFIter mfi(*averages[iavg], TileNoZ()); mfi.isValid(); ++mfi)
822  {
823  Box pbx = mfi.tilebox();
824  pbx.setSmall(2,0); pbx.setBig(2,0);
825 
826  const Array4<Real const>& T_mf_arr = fields[2]->const_array(mfi);
827  const Array4<Real const>& qv_mf_arr = (fields[3])? fields[3]->const_array(mfi) : Array4<const Real>{};
828  const Array4<Real const>& qr_mf_arr = (fields[4])? fields[4]->const_array(mfi) : Array4<const Real>{};
829 
830  if (m_interp) {
831  const auto plo = m_geom[lev].ProbLoArray();
832  const auto dxInv = m_geom[lev].InvCellSizeArray();
833  const auto z_phys_arr = z_phys->const_array(mfi);
834  auto x_pos_arr = x_pos->array(mfi);
835  auto y_pos_arr = y_pos->array(mfi);
836  auto z_pos_arr = z_pos->array(mfi);
837  ParallelFor(Gpu::KernelInfo().setReduction(true), pbx, [=]
838  AMREX_GPU_DEVICE(int i, int j, int k, Gpu::Handler const& handler) noexcept
839  {
840  Real T_interp{0};
841  Real qv_interp{0};
842  trilinear_interp_T(x_pos_arr(i,j,k), y_pos_arr(i,j,k), z_pos_arr(i,j,k),
843  &T_interp, T_mf_arr, z_phys_arr, plo, dxInv, 1);
844  trilinear_interp_T(x_pos_arr(i,j,k), y_pos_arr(i,j,k), z_pos_arr(i,j,k),
845  &qv_interp, qv_mf_arr, z_phys_arr, plo, dxInv, 1);
846  Real vfac;
847  if (qr_mf_arr) {
848  // We also have liquid water
849  Real qr_interp{0};
850  trilinear_interp_T(x_pos_arr(i,j,k), y_pos_arr(i,j,k), z_pos_arr(i,j,k),
851  &qr_interp, qr_mf_arr, z_phys_arr, plo, dxInv, 1);
852  vfac = 1.0 + 0.61*qv_interp - qr_interp;
853  } else {
854  vfac = 1.0 + 0.61*qv_interp;
855  }
856  const Real val = T_interp * vfac;
857  Gpu::deviceReduceSum(&plane_avg[iavg], val, handler);
858  });
859  } else {
860  auto k_arr = k_indx->const_array(mfi);
861  auto j_arr = j_indx ? j_indx->const_array(mfi) : Array4<const int> {};
862  auto i_arr = i_indx ? i_indx->const_array(mfi) : Array4<const int> {};
863  ParallelFor(Gpu::KernelInfo().setReduction(true), pbx, [=]
864  AMREX_GPU_DEVICE(int i, int j, int k, Gpu::Handler const& handler) noexcept
865  {
866  int mk = k_arr(i,j,k);
867  int mj = j_arr ? j_arr(i,j,k) : j;
868  int mi = i_arr ? i_arr(i,j,k) : i;
869  Real vfac;
870  if (qr_mf_arr) {
871  // We also have liquid water
872  vfac = 1.0 + 0.61*qv_mf_arr(mi,mj,mk) - qr_mf_arr(mi,mj,mk);
873  } else {
874  vfac = 1.0 + 0.61*qv_mf_arr(mi,mj,mk);
875  }
876  const Real val = T_mf_arr(mi,mj,mk) * vfac;
877  Gpu::deviceReduceSum(&plane_avg[iavg], val, handler);
878  });
879  }
880  }
881  }
882  else // copy temperature
883  {
884  int iavg = m_navg - 2;
885  denom[iavg] = 1.0 / (Real)ncell_plane[iavg];
886  // plane_avg[iavg] = plane_avg[2]
887  Gpu::copy(Gpu::deviceToDevice, pavg.begin() + 2, pavg.begin() + 3,
888  pavg.begin() + iavg);
889  }
890 
891  //
892  //------------------------------------------------------------------------
893  // Averages for the tangential velocity magnitude
894  // (This is cell-centered so we don't need to worry about double-counting)
895  //------------------------------------------------------------------------
896  //
897  {
898  int imf = 0;
899  int iavg = m_navg - 1;
900  denom[iavg] = 1.0 / (Real)ncell_plane[iavg];
901  val_old[iavg] = plane_average[iavg]*d_fact_old;
902 
903  const Real Vsg = m_Vsg[lev];
904 
905 #ifdef _OPENMP
906 #pragma omp parallel if (Gpu::notInLaunchRegion())
907 #endif
908  for (MFIter mfi(*averages[iavg], TileNoZ()); mfi.isValid(); ++mfi)
909  {
910  Box pbx = mfi.tilebox();
911  pbx.setSmall(2,0); pbx.setBig(2,0);
912 
913  // Last element is Umag and always cell centered
914  auto u_mf_arr = (m_rotate) ? rot_fields[imf ]->const_array(mfi) :
915  fields[imf ]->const_array(mfi);
916  auto v_mf_arr = (m_rotate) ? rot_fields[imf+1]->const_array(mfi) :
917  fields[imf+1]->const_array(mfi);
918 
919  if (m_interp) {
920  const auto plo = m_geom[lev].ProbLoArray();
921  const auto dxInv = m_geom[lev].InvCellSizeArray();
922  const auto z_phys_arr = z_phys->const_array(mfi);
923  auto x_pos_arr = x_pos->array(mfi);
924  auto y_pos_arr = y_pos->array(mfi);
925  auto z_pos_arr = z_pos->array(mfi);
926  ParallelFor(Gpu::KernelInfo().setReduction(true), pbx, [=]
927  AMREX_GPU_DEVICE(int i, int j, int k, Gpu::Handler const& handler) noexcept
928  {
929  Real u_interp{0};
930  Real v_interp{0};
931  trilinear_interp_T(x_pos_arr(i,j,k), y_pos_arr(i,j,k), z_pos_arr(i,j,k),
932  &u_interp, u_mf_arr, z_phys_arr, plo, dxInv, 1);
933  trilinear_interp_T(x_pos_arr(i,j,k), y_pos_arr(i,j,k), z_pos_arr(i,j,k),
934  &v_interp, v_mf_arr, z_phys_arr, plo, dxInv, 1);
935  const Real val = std::sqrt(u_interp*u_interp + v_interp*v_interp + Vsg*Vsg);
936  Gpu::deviceReduceSum(&plane_avg[iavg], val, handler);
937  });
938  } else {
939  auto k_arr = k_indx->const_array(mfi);
940  auto j_arr = j_indx ? j_indx->const_array(mfi) : Array4<const int> {};
941  auto i_arr = i_indx ? i_indx->const_array(mfi) : Array4<const int> {};
942  ParallelFor(Gpu::KernelInfo().setReduction(true), pbx, [=]
943  AMREX_GPU_DEVICE(int i, int j, int k, Gpu::Handler const& handler) noexcept
944  {
945  int mk = k_arr(i,j,k);
946  int mj = j_arr ? j_arr(i,j,k) : j;
947  int mi = i_arr ? i_arr(i,j,k) : i;
948  const Real u_val = 0.5 * (u_mf_arr(mi,mj,mk) + u_mf_arr(mi+1,mj ,mk));
949  const Real v_val = 0.5 * (v_mf_arr(mi,mj,mk) + v_mf_arr(mi ,mj+1,mk));
950  const Real val = std::sqrt(u_val*u_val + v_val*v_val + Vsg*Vsg);
951  Gpu::deviceReduceSum(&plane_avg[iavg], val, handler);
952  });
953  }
954  }
955  }
956 
957  // Copy to host and sum across procs
958  Gpu::copy(Gpu::deviceToHost, pavg.begin(), pavg.end(), plane_average.begin());
959  ParallelDescriptor::ReduceRealSum(plane_average.data(), plane_average.size());
960 
961  // No spatial variation with plane averages
962  for (int iavg(0); iavg < m_navg; ++iavg){
963  plane_average[iavg] *= denom[iavg]*d_fact_new;
964  plane_average[iavg] += val_old[iavg];
965  averages[iavg]->setVal(plane_average[iavg]);
966  }
967 }
TileNoZ
AMREX_FORCE_INLINE amrex::IntVect TileNoZ()
Definition: ERF_TileNoZ.H:11
MOSTAverage::m_navg
int m_navg
Definition: ERF_MOSTAverage.H:186
MOSTAverage::m_averages
amrex::Vector< amrex::Vector< std::unique_ptr< amrex::MultiFab > > > m_averages
Definition: ERF_MOSTAverage.H:197
MOSTAverage::m_y_pos
amrex::Vector< std::unique_ptr< amrex::MultiFab > > m_y_pos
Definition: ERF_MOSTAverage.H:192
MOSTAverage::m_i_indx
amrex::Vector< std::unique_ptr< amrex::iMultiFab > > m_i_indx
Definition: ERF_MOSTAverage.H:194
MOSTAverage::m_z_phys_nd
amrex::Vector< amrex::MultiFab * > m_z_phys_nd
Definition: ERF_MOSTAverage.H:179
MOSTAverage::m_x_pos
amrex::Vector< std::unique_ptr< amrex::MultiFab > > m_x_pos
Definition: ERF_MOSTAverage.H:191
MOSTAverage::m_Vsg
amrex::Vector< amrex::Real > m_Vsg
Definition: ERF_MOSTAverage.H:226
MOSTAverage::m_rot_fields
amrex::Vector< amrex::Vector< std::unique_ptr< amrex::MultiFab > > > m_rot_fields
Definition: ERF_MOSTAverage.H:198
MOSTAverage::m_plane_average
amrex::Vector< amrex::Vector< amrex::Real > > m_plane_average
Definition: ERF_MOSTAverage.H:203
MOSTAverage::m_z_pos
amrex::Vector< std::unique_ptr< amrex::MultiFab > > m_z_pos
Definition: ERF_MOSTAverage.H:193
MOSTAverage::m_ncell_plane
amrex::Vector< amrex::Vector< int > > m_ncell_plane
Definition: ERF_MOSTAverage.H:202
MOSTAverage::m_fact_new
amrex::Real m_fact_new
Definition: ERF_MOSTAverage.H:221
MOSTAverage::m_j_indx
amrex::Vector< std::unique_ptr< amrex::iMultiFab > > m_j_indx
Definition: ERF_MOSTAverage.H:195
MOSTAverage::m_fields
amrex::Vector< amrex::Vector< amrex::MultiFab * > > m_fields
Definition: ERF_MOSTAverage.H:178
MOSTAverage::m_fact_old
amrex::Real m_fact_old
Definition: ERF_MOSTAverage.H:221
MOSTAverage::m_interp
bool m_interp
Definition: ERF_MOSTAverage.H:213
MOSTAverage::m_geom
const amrex::Vector< amrex::Geometry > m_geom
Definition: ERF_MOSTAverage.H:177
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:118
MOSTAverage::m_k_indx
amrex::Vector< std::unique_ptr< amrex::iMultiFab > > m_k_indx
Definition: ERF_MOSTAverage.H:196

Referenced by compute_averages().

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

◆ compute_region_averages()

void MOSTAverage::compute_region_averages ( const int &  lev)

Function to compute average over local region.

Parameters
[in]levCurrent level
977 {
978  // Peel back the level
979  auto& fields = m_fields[lev];
980  auto& rot_fields = m_rot_fields[lev];
981  auto& averages = m_averages[lev];
982  const auto & geom = m_geom[lev];
983 
984  auto& z_phys = m_z_phys_nd[lev];
985  auto& x_pos = m_x_pos[lev];
986  auto& y_pos = m_y_pos[lev];
987  auto& z_pos = m_z_pos[lev];
988 
989  auto& i_indx = m_i_indx[lev];
990  auto& j_indx = m_j_indx[lev];
991  auto& k_indx = m_k_indx[lev];
992 
993  // Set factors for time averaging
994  Real d_fact_new, d_fact_old;
995  if (m_t_avg && m_t_init[lev]) {
996  d_fact_new = m_fact_new;
997  d_fact_old = m_fact_old;
998  } else {
999  d_fact_new = 1.0;
1000  d_fact_old = 0.0;
1001  }
1002 
1003  // Number of cells contained in the local average
1004  const Real denom = 1.0 / (Real) m_ncell_region;
1005 
1006  // Capture radius for device
1007  int d_radius = m_radius;
1008 
1009  //
1010  //----------------------------------------------------------
1011  // Averages for U,V,T,Qv
1012  //----------------------------------------------------------
1013  //
1014  for (int imf(0); imf < 4; ++imf) {
1015 
1016  // Continue if no valid Qv pointer
1017  if (!fields[imf]) continue;
1018 
1019 #ifdef _OPENMP
1020 #pragma omp parallel if (Gpu::notInLaunchRegion())
1021 #endif
1022  for (MFIter mfi(*fields[imf], TileNoZ()); mfi.isValid(); ++mfi) {
1023  Box pbx = mfi.tilebox(); pbx.setSmall(2,0); pbx.setBig(2,0);
1024 
1025  auto mf_arr = (m_rotate) ? rot_fields[imf]->const_array(mfi) :
1026  fields[imf]->const_array(mfi);
1027  auto ma_arr = averages[imf]->array(mfi);
1028 
1029  if (m_interp) {
1030  const auto plo = geom.ProbLoArray();
1031  const auto dx = geom.CellSizeArray();
1032  const auto dxInv = geom.InvCellSizeArray();
1033  const auto z_phys_arr = z_phys->const_array(mfi);
1034  auto x_pos_arr = x_pos->array(mfi);
1035  auto y_pos_arr = y_pos->array(mfi);
1036  auto z_pos_arr = z_pos->array(mfi);
1037  ParallelFor(pbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
1038  {
1039  ma_arr(i,j,k) *= d_fact_old;
1040 
1041  Real met_h_zeta = Compute_h_zeta_AtCellCenter(i,j,k,dxInv,z_phys_arr);
1042  for (int lk(-d_radius); lk <= (d_radius); ++lk) {
1043  for (int lj(-d_radius); lj <= (d_radius); ++lj) {
1044  for (int li(-d_radius); li <= (d_radius); ++li) {
1045  Real interp{0};
1046  Real xp = x_pos_arr(i+li,j+lj,k);
1047  Real yp = y_pos_arr(i+li,j+lj,k);
1048  Real zp = z_pos_arr(i+li,j+lj,k) + met_h_zeta*lk*dx[2];
1049  trilinear_interp_T(xp, yp, zp, &interp, mf_arr, z_phys_arr, plo, dxInv, 1);
1050  Real val = denom * interp * d_fact_new;
1051  ma_arr(i,j,k) += val;
1052  }
1053  }
1054  }
1055  });
1056  } else {
1057  auto k_arr = k_indx->const_array(mfi);
1058  auto j_arr = j_indx ? j_indx->const_array(mfi) : Array4<const int> {};
1059  auto i_arr = i_indx ? i_indx->const_array(mfi) : Array4<const int> {};
1060  ParallelFor(pbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
1061  {
1062  ma_arr(i,j,k) *= d_fact_old;
1063 
1064  int mk = k_arr(i,j,k);
1065  int mj = j_arr ? j_arr(i,j,k) : j;
1066  int mi = i_arr ? i_arr(i,j,k) : i;
1067  for (int lk(mk-d_radius); lk <= (mk+d_radius); ++lk) {
1068  for (int lj(mj-d_radius); lj <= (mj+d_radius); ++lj) {
1069  for (int li(mi-d_radius); li <= (mi+d_radius); ++li) {
1070  Real val = denom * mf_arr(li, lj, lk) * d_fact_new;
1071  ma_arr(i,j,k) += val;
1072  }
1073  }
1074  }
1075  });
1076  }
1077  } // MFiter
1078 
1079  // Fill interior ghost cells and any ghost cells outside a periodic domain
1080  //***********************************************************************************
1081  averages[imf]->FillBoundary(geom.periodicity());
1082 
1083  } // imf
1084 
1085  //
1086  //----------------------------------------------------------
1087  // Averages for virtual potential temperature
1088  //----------------------------------------------------------
1089  //
1090  if (fields[3]) // We have water vapor
1091  {
1092  int iavg = 4;
1093 
1094 #ifdef _OPENMP
1095 #pragma omp parallel if (Gpu::notInLaunchRegion())
1096 #endif
1097  for (MFIter mfi(*averages[iavg], TileNoZ()); mfi.isValid(); ++mfi) {
1098  Box pbx = mfi.tilebox(); pbx.setSmall(2,0); pbx.setBig(2,0);
1099 
1100  const Array4<Real const>& T_mf_arr = fields[2]->const_array(mfi);
1101  const Array4<Real const>& qv_mf_arr = (fields[3])? fields[3]->const_array(mfi) : Array4<const Real>{};
1102  const Array4<Real const>& qr_mf_arr = (fields[4])? fields[4]->const_array(mfi) : Array4<const Real>{};
1103  auto ma_arr = averages[iavg]->array(mfi);
1104 
1105  if (m_interp) {
1106  const auto plo = geom.ProbLoArray();
1107  const auto dx = geom.CellSizeArray();
1108  const auto dxInv = geom.InvCellSizeArray();
1109  const auto z_phys_arr = z_phys->const_array(mfi);
1110  auto x_pos_arr = x_pos->array(mfi);
1111  auto y_pos_arr = y_pos->array(mfi);
1112  auto z_pos_arr = z_pos->array(mfi);
1113  ParallelFor(pbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
1114  {
1115  ma_arr(i,j,k) *= d_fact_old;
1116 
1117  Real met_h_zeta = Compute_h_zeta_AtCellCenter(i,j,k,dxInv,z_phys_arr);
1118  for (int lk(-d_radius); lk <= (d_radius); ++lk) {
1119  for (int lj(-d_radius); lj <= (d_radius); ++lj) {
1120  for (int li(-d_radius); li <= (d_radius); ++li) {
1121  Real T_interp{0};
1122  Real qv_interp{0};
1123  Real xp = x_pos_arr(i+li,j+lj,k);
1124  Real yp = y_pos_arr(i+li,j+lj,k);
1125  Real zp = z_pos_arr(i+li,j+lj,k) + met_h_zeta*lk*dx[2];
1126  trilinear_interp_T(xp, yp, zp, &T_interp, T_mf_arr, z_phys_arr, plo, dxInv, 1);
1127  trilinear_interp_T(xp, yp, zp, &qv_interp, qv_mf_arr, z_phys_arr, plo, dxInv, 1);
1128  Real vfac;
1129  if (qr_mf_arr) {
1130  // We also have liquid water
1131  Real qr_interp{0};
1132  trilinear_interp_T(x_pos_arr(i,j,k), y_pos_arr(i,j,k), z_pos_arr(i,j,k),
1133  &qr_interp, qr_mf_arr, z_phys_arr, plo, dxInv, 1);
1134  vfac = 1.0 + 0.61*qv_interp - qr_interp;
1135  } else {
1136  vfac = 1.0 + 0.61*qv_interp;
1137  }
1138  const Real mag = T_interp * vfac;
1139  const Real val = denom * mag * d_fact_new;
1140  ma_arr(i,j,k) += val;
1141  }
1142  }
1143  }
1144  });
1145  } else {
1146  auto k_arr = k_indx->const_array(mfi);
1147  auto j_arr = j_indx ? j_indx->const_array(mfi) : Array4<const int> {};
1148  auto i_arr = i_indx ? i_indx->const_array(mfi) : Array4<const int> {};
1149  ParallelFor(pbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
1150  {
1151  ma_arr(i,j,k) *= d_fact_old;
1152 
1153  int mk = k_arr(i,j,k);
1154  int mj = j_arr ? j_arr(i,j,k) : j;
1155  int mi = i_arr ? i_arr(i,j,k) : i;
1156  for (int lk(mk-d_radius); lk <= (mk+d_radius); ++lk) {
1157  for (int lj(mj-d_radius); lj <= (mj+d_radius); ++lj) {
1158  for (int li(mi-d_radius); li <= (mi+d_radius); ++li) {
1159  Real vfac;
1160  if (qr_mf_arr) {
1161  // We also have liquid water
1162  vfac = 1.0 + 0.61*qv_mf_arr(li,lj,lk) - qr_mf_arr(li,lj,lk);
1163  } else {
1164  vfac = 1.0 + 0.61*qv_mf_arr(li,lj,lk);
1165  }
1166  const Real mag = T_mf_arr(li,lj,lk) * vfac;
1167  const Real val = denom * mag * d_fact_new;
1168  ma_arr(i,j,k) += val;
1169  }
1170  }
1171  }
1172  });
1173  }
1174  } // MFiter
1175 
1176  // Fill interior ghost cells and any ghost cells outside a periodic domain
1177  //***********************************************************************************
1178  averages[iavg]->FillBoundary(geom.periodicity());
1179 
1180  }
1181  else // copy temperature
1182  {
1183  int iavg = m_navg - 2;
1184  IntVect ng = averages[iavg]->nGrowVect();
1185  MultiFab::Copy(*(averages[iavg]),*(averages[2]),0,0,1,ng);
1186  }
1187 
1188  //
1189  //----------------------------------------------------------
1190  // Averages for the tangential velocity magnitude
1191  //----------------------------------------------------------
1192  //
1193  {
1194  int imf = 0;
1195  int iavg = m_navg - 1;
1196 
1197  const Real Vsg = m_Vsg[lev];
1198 
1199 #ifdef _OPENMP
1200 #pragma omp parallel if (Gpu::notInLaunchRegion())
1201 #endif
1202  for (MFIter mfi(*averages[iavg], TileNoZ()); mfi.isValid(); ++mfi) {
1203  Box pbx = mfi.tilebox(); pbx.setSmall(2,0); pbx.setBig(2,0);
1204 
1205  auto u_mf_arr = (m_rotate) ? rot_fields[imf ]->const_array(mfi) :
1206  fields[imf ]->const_array(mfi);
1207  auto v_mf_arr = (m_rotate) ? rot_fields[imf+1]->const_array(mfi) :
1208  fields[imf+1]->const_array(mfi);
1209  auto ma_arr = averages[iavg]->array(mfi);
1210 
1211  if (m_interp) {
1212  const auto plo = geom.ProbLoArray();
1213  const auto dx = geom.CellSizeArray();
1214  const auto dxInv = geom.InvCellSizeArray();
1215  const auto z_phys_arr = z_phys->const_array(mfi);
1216  auto x_pos_arr = x_pos->array(mfi);
1217  auto y_pos_arr = y_pos->array(mfi);
1218  auto z_pos_arr = z_pos->array(mfi);
1219  ParallelFor(pbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
1220  {
1221  ma_arr(i,j,k) *= d_fact_old;
1222 
1223  Real met_h_zeta = Compute_h_zeta_AtCellCenter(i,j,k,dxInv,z_phys_arr);
1224  for (int lk(-d_radius); lk <= (d_radius); ++lk) {
1225  for (int lj(-d_radius); lj <= (d_radius); ++lj) {
1226  for (int li(-d_radius); li <= (d_radius); ++li) {
1227  Real u_interp{0};
1228  Real v_interp{0};
1229  Real xp = x_pos_arr(i+li,j+lj,k);
1230  Real yp = y_pos_arr(i+li,j+lj,k);
1231  Real zp = z_pos_arr(i+li,j+lj,k) + met_h_zeta*lk*dx[2];
1232  trilinear_interp_T(xp, yp, zp, &u_interp, u_mf_arr, z_phys_arr, plo, dxInv, 1);
1233  trilinear_interp_T(xp, yp, zp, &v_interp, v_mf_arr, z_phys_arr, plo, dxInv, 1);
1234  const Real mag = std::sqrt(u_interp*u_interp + v_interp*v_interp + Vsg*Vsg);
1235  Real val = denom * mag * d_fact_new;
1236  ma_arr(i,j,k) += val;
1237  }
1238  }
1239  }
1240  });
1241  } else {
1242  auto k_arr = k_indx->const_array(mfi);
1243  auto j_arr = j_indx ? j_indx->const_array(mfi) : Array4<const int> {};
1244  auto i_arr = i_indx ? i_indx->const_array(mfi) : Array4<const int> {};
1245  ParallelFor(pbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
1246  {
1247  ma_arr(i,j,k) *= d_fact_old;
1248 
1249  int mk = k_arr(i,j,k);
1250  int mj = j_arr ? j_arr(i,j,k) : j;
1251  int mi = i_arr ? i_arr(i,j,k) : i;
1252  for (int lk(mk-d_radius); lk <= (mk+d_radius); ++lk) {
1253  for (int lj(mj-d_radius); lj <= (mj+d_radius); ++lj) {
1254  for (int li(mi-d_radius); li <= (mi+d_radius); ++li) {
1255  const Real u_val = 0.5 * (u_mf_arr(li,lj,lk) + u_mf_arr(li+1,lj ,lk));
1256  const Real v_val = 0.5 * (v_mf_arr(li,lj,lk) + v_mf_arr(li ,lj+1,lk));
1257  const Real mag = std::sqrt(u_val*u_val + v_val*v_val + Vsg*Vsg);
1258  Real val = denom * mag * d_fact_new;
1259  ma_arr(i,j,k) += val;
1260  }
1261  }
1262  }
1263  });
1264  }
1265  } // MFiter
1266 
1267  // Fill interior ghost cells and any ghost cells outside a periodic domain
1268  //***********************************************************************************
1269  averages[iavg]->FillBoundary(geom.periodicity());
1270 
1271  }
1272 
1273  // NOTE: Checking periodicity with the geom structure is not
1274  // sufficient at higher levels. The BA may be contained
1275  // within the domain and it's exterior ghost cells filled
1276  // from interpolation; yet the domain BCs are periodic.
1277 
1278  // Need to fill ghost cells outside the domain if not periodic
1279  bool not_per_x = !(geom.periodicity().isPeriodic(0));
1280  bool not_per_y = !(geom.periodicity().isPeriodic(1));
1281  Box cc_bnd_bx = (m_fields[lev][2]->boxArray()).minimalBox();
1282  Box domain = geom.Domain();
1283  if (domain.contains(cc_bnd_bx) || (not_per_x || not_per_y)) {
1284  for (int iavg(0); iavg < m_navg; ++iavg) {
1285  IntVect ng = averages[iavg]->nGrowVect(); ng[2]=0;
1286 
1287  // NOTE: Level 0 spans the whole domain, but finer
1288  // levels do not have such a restriction.
1289  // For now, use the bounding box of the boxArray.
1290 
1291  // NOTE2: The fields and averages have different indexing.
1292  // The averages are: U/V/T/Qv/Tv/Umag
1293  // The fields are: U/V/T/Qv/Qr/W
1294  // We clip iavg at 2 since all the remaining data is CC
1295 
1296  // Bounded box of CC data used for normalization
1297  int imf = min(iavg,2);
1298  Box bnd_bx = (m_fields[lev][imf]->boxArray()).minimalBox();
1299 #ifdef _OPENMP
1300 #pragma omp parallel if (Gpu::notInLaunchRegion())
1301 #endif
1302  for (MFIter mfi(*averages[iavg], TileNoZ()); mfi.isValid(); ++mfi) {
1303  Box gpbx = mfi.growntilebox(ng); gpbx.setSmall(2,0); gpbx.setBig(2,0);
1304 
1305  if (bnd_bx.contains(gpbx)) continue;
1306 
1307  auto ma_arr = averages[iavg]->array(mfi);
1308 
1309  int i_lo = bnd_bx.smallEnd(0); int i_hi = bnd_bx.bigEnd(0);
1310  int j_lo = bnd_bx.smallEnd(1); int j_hi = bnd_bx.bigEnd(1);
1311  ParallelFor(gpbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
1312  {
1313  int li, lj;
1314  li = i < i_lo ? i_lo : i;
1315  li = li > i_hi ? i_hi : li;
1316  lj = j < j_lo ? j_lo : j;
1317  lj = lj > j_hi ? j_hi : lj;
1318 
1319  ma_arr(i,j,k) = ma_arr(li,lj,k);
1320  });
1321  } // MFiter
1322  } // iavg
1323  } // Not periodic
1324 }
Compute_h_zeta_AtCellCenter
AMREX_FORCE_INLINE AMREX_GPU_DEVICE amrex::Real Compute_h_zeta_AtCellCenter(const int &i, const int &j, const int &k, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &z_nd)
Definition: ERF_TerrainMetrics.H:46
MOSTAverage::m_radius
int m_radius
Definition: ERF_MOSTAverage.H:207
MOSTAverage::m_ncell_region
int m_ncell_region
Definition: ERF_MOSTAverage.H:208

Referenced by compute_averages().

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

◆ get_average()

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

Referenced by ABLMost::get_mac_avg().

Here is the caller graph for this function:

◆ get_zref()

amrex::Real MOSTAverage::get_zref ( ) const
inline
103 { return m_zref; }
MOSTAverage::m_zref
amrex::Real m_zref
Definition: ERF_MOSTAverage.H:190

Referenced by ABLMost::get_zref().

Here is the caller graph for this function:

◆ make_MOSTAverage_at_level()

void MOSTAverage::make_MOSTAverage_at_level ( const int &  lev,
const amrex::Vector< amrex::MultiFab * > &  vars_old,
std::unique_ptr< amrex::MultiFab > &  Theta_prim,
std::unique_ptr< amrex::MultiFab > &  Qv_prim,
std::unique_ptr< amrex::MultiFab > &  Qr_prim,
std::unique_ptr< amrex::MultiFab > &  z_phys_nd 
)
74 {
75  m_fields[lev].resize(m_nvar);
76  m_rot_fields[lev].resize(m_nvar-1);
77  m_averages[lev].resize(m_navg);
78  m_z_phys_nd[lev] = z_phys_nd.get();
79 
80  bool use_terrain_fitted_coords = ( (m_terrain_type == TerrainType::StaticFittedMesh) ||
81  (m_terrain_type == TerrainType::MovingFittedMesh) );
82 
83  { // Nodal in x
84  auto& mf = *vars_old[Vars::xvel];
85  // Create a 2D ba, dm, & ghost cells
86  const BoxArray& ba = mf.boxArray();
87  BoxList bl2d = ba.boxList();
88  for (auto& b : bl2d) b.setRange(2,0);
89  BoxArray ba2d(std::move(bl2d));
90  const DistributionMapping& dm = mf.DistributionMap();
91  const int ncomp = 1;
92  IntVect ng = mf.nGrowVect(); ng[2]=0;
93 
94  m_fields[lev][0] = vars_old[Vars::xvel];
95  m_averages[lev][0] = std::make_unique<MultiFab>(ba2d,dm,ncomp,ng);
96  m_averages[lev][0]->setVal(1.E34);
97  if (m_rotate) {
98  m_rot_fields[lev][0] = std::make_unique<MultiFab>(ba,dm,ncomp,ng);
99  MultiFab::Copy(*m_rot_fields[lev][0],mf,0,0,1,ng);
100  } else {
101  m_rot_fields[lev][0] = nullptr;
102  }
103  }
104  { // Nodal in y
105  auto& mf = *vars_old[Vars::yvel];
106  // Create a 2D ba, dm, & ghost cells
107  const BoxArray& ba = mf.boxArray();
108  BoxList bl2d = ba.boxList();
109  for (auto& b : bl2d) b.setRange(2,0);
110  BoxArray ba2d(std::move(bl2d));
111  const DistributionMapping& dm = mf.DistributionMap();
112  const int ncomp = 1;
113  IntVect ng = mf.nGrowVect(); ng[2]=0;
114 
115  m_fields[lev][1] = vars_old[Vars::yvel];
116  m_averages[lev][1] = std::make_unique<MultiFab>(ba2d,dm,ncomp,ng);
117  m_averages[lev][1]->setVal(1.E34);
118  if (m_rotate) {
119  m_rot_fields[lev][1] = std::make_unique<MultiFab>(ba,dm,ncomp,ng);
120  MultiFab::Copy(*m_rot_fields[lev][1],mf,0,0,1,ng);
121  } else {
122  m_rot_fields[lev][1] = nullptr;
123  }
124  }
125  { // CC vars
126  auto& mf = *Theta_prim;
127  // Create a 2D ba, dm, & ghost cells
128  const BoxArray& ba = mf.boxArray();
129  BoxList bl2d = ba.boxList();
130  for (auto& b : bl2d) b.setRange(2,0);
131  BoxArray ba2d(std::move(bl2d));
132  const DistributionMapping& dm = mf.DistributionMap();
133  const int ncomp = 1;
134  const int incomp = 1;
135  IntVect ng = mf.nGrowVect(); ng[2]=0;
136 
137  // Get field pointers
138  m_fields[lev][2] = Theta_prim.get();
139  m_fields[lev][3] = Qv_prim.get();
140  m_fields[lev][4] = Qr_prim.get();
141 
142  // Initialize remaining multifabs
143  for (int iavg(2); iavg < m_navg; ++iavg) {
144  m_averages[lev][iavg] = std::make_unique<MultiFab>(ba2d,dm,ncomp,ng);
145  m_averages[lev][iavg]->setVal(1.E34);
146  }
147 
148  // Default to dry
149  m_averages[lev][3]->setVal(0.0);
150 
151  if (m_rotate) {
152  m_rot_fields[lev][2] = std::make_unique<MultiFab>(ba,dm,ncomp,ng);
153  m_rot_fields[lev][3] = std::make_unique<MultiFab>(ba,dm,ncomp,ng);
154  MultiFab::Copy(*m_rot_fields[lev][2],*Theta_prim,0,0,1,ng);
155  if (Qv_prim) MultiFab::Copy(*m_rot_fields[lev][3],*Qv_prim,0,0,1,ng);
156  } else {
157  m_rot_fields[lev][2] = nullptr;
158  m_rot_fields[lev][3] = nullptr;
159  }
160 
161  if (use_terrain_fitted_coords && m_norm_vec && m_interp) {
162  m_x_pos[lev] = std::make_unique<MultiFab>(ba2d,dm,ncomp,ng);
163  m_y_pos[lev] = std::make_unique<MultiFab>(ba2d,dm,ncomp,ng);
164  m_z_pos[lev] = std::make_unique<MultiFab>(ba2d,dm,ncomp,ng);
165  } else if (use_terrain_fitted_coords && m_interp) {
166  m_x_pos[lev] = std::make_unique<MultiFab>(ba2d,dm,ncomp,ng);
167  m_y_pos[lev] = std::make_unique<MultiFab>(ba2d,dm,ncomp,ng);
168  m_z_pos[lev] = std::make_unique<MultiFab>(ba2d,dm,ncomp,ng);
169  } else if (use_terrain_fitted_coords && m_norm_vec) {
170  m_i_indx[lev] = std::make_unique<iMultiFab>(ba2d,dm,incomp,ng);
171  m_j_indx[lev] = std::make_unique<iMultiFab>(ba2d,dm,incomp,ng);
172  m_k_indx[lev] = std::make_unique<iMultiFab>(ba2d,dm,incomp,ng);
173  } else {
174  m_k_indx[lev] = std::make_unique<iMultiFab>(ba2d,dm,incomp,ng);
175  }
176  }
177  // Nodal in z (only used with terrain stress rotations)
178  m_fields[lev][4] = vars_old[Vars::zvel];
179 
180  // Setup auxiliary data for spatial configuration & policy
181  //--------------------------------------------------------
182  if (use_terrain_fitted_coords && m_norm_vec && m_interp) { // Terrain w/ norm & w/ interpolation
183  set_norm_positions_T(lev);
184  } else if (use_terrain_fitted_coords && m_interp) { // Terrain w/ interpolation
185  set_z_positions_T(lev);
186  } else if (use_terrain_fitted_coords && m_norm_vec) { // Terrain w/ norm & w/o interpolation
187  set_norm_indices_T(lev);
188  } else if (use_terrain_fitted_coords) { // Terrain
189  set_k_indices_T(lev);
190  } else { // No Terrain
191  set_k_indices_N(lev);
192  }
193 
194  // Setup normalization data for the chosen policy
195  //--------------------------------------------------------
196  switch(m_policy) {
197  case 0: // Plane average
198  set_plane_normalization(lev);
199  break;
200  case 1: // Local region/point
201  set_region_normalization(lev);
202  break;
203  default:
204  AMREX_ASSERT_WITH_MESSAGE(false, "Unknown policy for MOSTAverage!");
205  }
206 
207  // Set up the exponential time filtering
208  //--------------------------------------------------------
209  if (m_t_avg) {
210  // Exponential filter function
211  m_fact_old = std::exp(-1.0 / m_time_window);
212 
213  // Enforce discrete normalization: (mfn*val_new + mfo*val_old)
214  m_fact_new = 1.0 - m_fact_old;
215 
216  // None of the averages are initialized
217  m_t_init.resize(m_maxlev,0);
218  }
219 
220  // Corrections to the mean surface velocity
221  m_Vsg = Vector<Real>(m_maxlev, 0.0);
222  if (include_subgrid_vel) {
223  if (include_subgrid_vel) {
224  Print() << "Subgrid velocity scale correction at level : " << lev << ' ';
225  const auto dxArr = m_geom[lev].CellSizeArray();
226  Real dx = std::sqrt(dxArr[0]*dxArr[1]);
227  if (dx > 5000.) {
228  m_Vsg[lev] = 0.32 * std::pow(dx/5000.-1, 0.33);
229  }
230  Print() << m_Vsg[lev] << std::endl;
231  }
232  }
233 }
MOSTAverage::set_region_normalization
void set_region_normalization(const int &)
Definition: ERF_MOSTAverage.H:59
MOSTAverage::set_z_positions_T
void set_z_positions_T(const int &lev)
Definition: ERF_MOSTAverage.cpp:561
MOSTAverage::set_norm_positions_T
void set_norm_positions_T(const int &lev)
Definition: ERF_MOSTAverage.cpp:605
MOSTAverage::set_k_indices_T
void set_k_indices_T(const int &lev)
Definition: ERF_MOSTAverage.cpp:435
MOSTAverage::set_plane_normalization
void set_plane_normalization(const int &lev)
Definition: ERF_MOSTAverage.cpp:322
MOSTAverage::m_terrain_type
TerrainType m_terrain_type
Definition: ERF_MOSTAverage.H:181
MOSTAverage::m_norm_vec
bool m_norm_vec
Definition: ERF_MOSTAverage.H:214
MOSTAverage::m_nvar
int m_nvar
Definition: ERF_MOSTAverage.H:185
MOSTAverage::set_norm_indices_T
void set_norm_indices_T(const int &lev)
Definition: ERF_MOSTAverage.cpp:490
MOSTAverage::m_maxlev
int m_maxlev
Definition: ERF_MOSTAverage.H:187
MOSTAverage::include_subgrid_vel
bool include_subgrid_vel
Definition: ERF_MOSTAverage.H:225
MOSTAverage::m_time_window
amrex::Real m_time_window
Definition: ERF_MOSTAverage.H:220
MOSTAverage::set_k_indices_N
void set_k_indices_N(const int &lev)
Definition: ERF_MOSTAverage.cpp:380
Vars::xvel
@ xvel
Definition: ERF_IndexDefines.H:141
Vars::zvel
@ zvel
Definition: ERF_IndexDefines.H:143
Vars::yvel
@ yvel
Definition: ERF_IndexDefines.H:142

Referenced by ABLMost::make_MOST_at_level().

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

◆ operator=() [1/2]

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

◆ operator=() [2/2]

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

◆ set_k_indices_N()

void MOSTAverage::set_k_indices_N ( const int &  lev)

Function to set K indices without terrain.

381 {
382  ParmParse pp(m_pp_prefix);
383  auto read_z = pp.query("most.zref",m_zref);
384  auto read_k = pp.queryarr("most.k_arr_in",m_k_in);
385 
386  // Default behavior is to use the first cell center
387  if (!read_z && !read_k) {
388  Real m_zlo = m_geom[0].ProbLo(2);
389  Real m_dz = m_geom[0].CellSize(2);
390  m_zref = m_zlo + 0.5 * m_dz;
391  Print() << "Reference height for MOST set to " << m_zref << std::endl;
392  read_z = true;
393  }
394 
395  // Specify z_ref & compute k_indx (z_ref takes precedence)
396  if (read_z) {
397  Real m_zlo = m_geom[lev].ProbLo(2);
398  Real m_zhi = m_geom[lev].ProbHi(2);
399  Real m_dz = m_geom[lev].CellSize(2);
400 
401  amrex::ignore_unused(m_zhi);
402 
403  AMREX_ASSERT_WITH_MESSAGE(m_zref >= m_zlo + 0.5 * m_dz,
404  "Query point must be past first z-cell!");
405 
406  AMREX_ASSERT_WITH_MESSAGE(m_zref <= m_zhi - 0.5 * m_dz,
407  "Query point must be below the last z-cell!");
408 
409  int lk = static_cast<int>(floor((m_zref - m_zlo) / m_dz - 0.5));
410 
411  m_zref = (lk + 0.5) * m_dz + m_zlo;
412 
413  AMREX_ALWAYS_ASSERT(lk >= m_radius);
414 
415  m_k_indx[lev]->setVal(lk);
416  // Specified k_indx & compute z_ref
417  } else if (read_k) {
418  AMREX_ASSERT_WITH_MESSAGE(m_k_in[lev] >= m_radius,
419  "K index must be larger than averaging radius!");
420  m_k_indx[lev]->setVal(m_k_in[lev]);
421 
422  // TODO: check that z_ref is constant across levels
423  Real m_zlo = m_geom[0].ProbLo(2);
424  Real m_dz = m_geom[0].CellSize(2);
425  m_zref = ((Real)m_k_in[0] + 0.5) * m_dz + m_zlo;
426  }
427 }
pp
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real pp(amrex::Real y)
Definition: ERF_MicrophysicsUtils.H:219
MOSTAverage::m_pp_prefix
std::string m_pp_prefix
Definition: ERF_MOSTAverage.H:180
MOSTAverage::m_k_in
amrex::Vector< int > m_k_in
Definition: ERF_MOSTAverage.H:209

Referenced by make_MOSTAverage_at_level().

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

◆ set_k_indices_T()

void MOSTAverage::set_k_indices_T ( const int &  lev)

Function to set K indices with terrain.

436 {
437  ParmParse pp(m_pp_prefix);
438  auto read_z = pp.query("most.zref",m_zref);
439  auto read_k = pp.queryarr("most.k_arr_in",m_k_in);
440 
441  // No default behavior with terrain (we can't tell the difference between
442  // vertical grid stretching and true terrain)
443  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(read_z != read_k,
444  "Need to specify zref or k_arr_in for MOST");
445 
446  // Capture for device
447  Real d_zref = m_zref;
448  Real d_radius = m_radius;
449  amrex::ignore_unused(d_radius);
450 
451  // Specify z_ref & compute k_indx (z_ref takes precedence)
452  if (read_z) {
453  int kmax = m_geom[lev].Domain().bigEnd(2);
454  for (MFIter mfi(*m_k_indx[lev], TileNoZ()); mfi.isValid(); ++mfi) {
455  Box npbx = mfi.tilebox(IntVect(1,1,0),IntVect(1,1,0));
456  const auto z_phys_arr = m_z_phys_nd[lev]->const_array(mfi);
457  auto k_arr = m_k_indx[lev]->array(mfi);
458  ParallelFor(npbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
459  {
460  k_arr(i,j,k) = 0;
461  Real z_bot_face = 0.25 * ( z_phys_arr(i ,j ,k) + z_phys_arr(i+1,j ,k)
462  + z_phys_arr(i ,j+1,k) + z_phys_arr(i+1,j+1,k) );
463  Real z_target = z_bot_face + d_zref;
464  for (int lk(0); lk<=kmax; ++lk) {
465  Real z_lo = 0.25 * ( z_phys_arr(i,j ,lk ) + z_phys_arr(i+1,j ,lk )
466  + z_phys_arr(i,j+1,lk ) + z_phys_arr(i+1,j+1,lk ) );
467  Real z_hi = 0.25 * ( z_phys_arr(i,j ,lk+1) + z_phys_arr(i+1,j ,lk+1)
468  + z_phys_arr(i,j+1,lk+1) + z_phys_arr(i+1,j+1,lk+1) );
469  if (z_target > z_lo && z_target < z_hi){
470  AMREX_ASSERT_WITH_MESSAGE(lk >= d_radius,
471  "K index must be larger than averaging radius!");
472  k_arr(i,j,k) = lk;
473  break;
474  }
475  }
476  });
477  }
478  // Specified k_indx & compute z_ref
479  } else if (read_k) {
480  AMREX_ASSERT_WITH_MESSAGE(false, "Specified k-indx with terrain not implemented!");
481  }
482 }

Referenced by make_MOSTAverage_at_level().

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

◆ set_norm_indices_T()

void MOSTAverage::set_norm_indices_T ( const int &  lev)

Function to set I,J,K indices with terrain normals.

491 {
492  ParmParse pp(m_pp_prefix);
493  pp.get("most.zref",m_zref);
494 
495  // Capture for device
496  Real d_zref = m_zref;
497  Real d_radius = m_radius;
498 
499  int kmax = m_geom[lev].Domain().bigEnd(2);
500  const auto dxInv = m_geom[lev].InvCellSizeArray();
501  IntVect ng = m_k_indx[lev]->nGrowVect(); ng[2]=0;
502  for (MFIter mfi(*m_k_indx[lev], TileNoZ()); mfi.isValid(); ++mfi) {
503  Box npbx = mfi.tilebox(IntVect(1,1,0),IntVect(1,1,0));
504  Box gpbx = mfi.growntilebox(ng);
505  const auto z_phys_arr = m_z_phys_nd[lev]->const_array(mfi);
506  auto i_arr = m_i_indx[lev]->array(mfi);
507  auto j_arr = m_j_indx[lev]->array(mfi);
508  auto k_arr = m_k_indx[lev]->array(mfi);
509  ParallelFor(npbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
510  {
511  // Elements of normal vector
512  Real met_h_xi = Compute_h_xi_AtCellCenter (i,j,k,dxInv,z_phys_arr);
513  Real met_h_eta = Compute_h_eta_AtCellCenter(i,j,k,dxInv,z_phys_arr);
514  Real mag = std::sqrt(met_h_xi*met_h_xi + met_h_eta*met_h_eta + 1.0);
515 
516  // Unit-normal vector scaled by z_ref
517  Real delta_x = -met_h_xi/mag * d_zref;
518  Real delta_y = -met_h_eta/mag * d_zref;
519  Real delta_z = 1.0/mag * d_zref;
520 
521  // Compute i & j as displacements (no grid stretching)
522  int delta_i = static_cast<int>(std::round(delta_x*dxInv[0]));
523  int delta_j = static_cast<int>(std::round(delta_y*dxInv[1]));
524  int i_new = i + delta_i;
525  int j_new = j + delta_j;
526  i_arr(i,j,k) = i_new;
527  j_arr(i,j,k) = j_new;
528 
529  // Search for k (grid is stretched in z)
530  Real z_bot_face = 0.25 * ( z_phys_arr(i ,j ,k) + z_phys_arr(i+1,j ,k)
531  + z_phys_arr(i ,j+1,k) + z_phys_arr(i+1,j+1,k) );
532  Real z_target = z_bot_face + delta_z;
533  for (int lk(0); lk<=kmax; ++lk) {
534  Real z_lo = 0.25 * ( z_phys_arr(i_new,j_new ,lk ) + z_phys_arr(i_new+1,j_new ,lk )
535  + z_phys_arr(i_new,j_new+1,lk ) + z_phys_arr(i_new+1,j_new+1,lk ) );
536  Real z_hi = 0.25 * ( z_phys_arr(i_new,j_new ,lk+1) + z_phys_arr(i_new+1,j_new ,lk+1)
537  + z_phys_arr(i_new,j_new+1,lk+1) + z_phys_arr(i_new+1,j_new+1,lk+1) );
538  if (z_target > z_lo && z_target < z_hi){
539  AMREX_ASSERT_WITH_MESSAGE(lk >= d_radius,
540  "K index must be larger than averaging radius!");
541  amrex::ignore_unused(d_radius);
542  k_arr(i,j,k) = lk;
543  break;
544  }
545  }
546 
547  // Destination cell must be contained on the current process!
548  amrex::ignore_unused(gpbx);
549  AMREX_ASSERT_WITH_MESSAGE(gpbx.contains(i_arr(i,j,k),j_arr(i,j,k),k_arr(i,j,k)),
550  "Query index outside of proc domain!");
551  });
552  }
553 }
Compute_h_eta_AtCellCenter
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_eta_AtCellCenter(const int &i, const int &j, const int &k, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &z_nd)
Definition: ERF_TerrainMetrics.H:76
Compute_h_xi_AtCellCenter
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_xi_AtCellCenter(const int &i, const int &j, const int &k, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &z_nd)
Definition: ERF_TerrainMetrics.H:61

Referenced by make_MOSTAverage_at_level().

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

◆ set_norm_positions_T()

void MOSTAverage::set_norm_positions_T ( const int &  lev)

Function to set positions with terrain and normal vector.

606 {
607  ParmParse pp(m_pp_prefix);
608  pp.get("most.zref",m_zref);
609 
610  // Capture for device
611  Real d_zref = m_zref;
612 
613  RealVect base;
614  const auto dx = m_geom[lev].CellSizeArray();
615  const auto dxInv = m_geom[lev].InvCellSizeArray();
616  IntVect ng = m_x_pos[lev]->nGrowVect(); ng[2]=0;
617  for (MFIter mfi(*m_x_pos[lev], TileNoZ()); mfi.isValid(); ++mfi) {
618  Box npbx = mfi.tilebox(IntVect(1,1,0),IntVect(1,1,0));
619  Box gpbx = mfi.growntilebox(ng);
620  RealBox grb{gpbx,dx.data(),base.dataPtr()};
621 
622  const auto z_phys_arr = m_z_phys_nd[lev]->const_array(mfi);
623  auto x_pos_arr = m_x_pos[lev]->array(mfi);
624  auto y_pos_arr = m_y_pos[lev]->array(mfi);
625  auto z_pos_arr = m_z_pos[lev]->array(mfi);
626  ParallelFor(npbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
627  {
628  // Elements of normal vector
629  Real met_h_xi = Compute_h_xi_AtCellCenter (i,j,k,dxInv,z_phys_arr);
630  Real met_h_eta = Compute_h_eta_AtCellCenter(i,j,k,dxInv,z_phys_arr);
631  Real mag = std::sqrt(met_h_xi*met_h_xi + met_h_eta*met_h_eta + 1.0);
632 
633  // Unit-normal vector scaled by z_ref
634  Real delta_x = -met_h_xi/mag * d_zref;
635  Real delta_y = -met_h_eta/mag * d_zref;
636  Real delta_z = 1.0/mag * d_zref;
637 
638  // Position of the current node (indx:0,0,1)
639  Real x0 = ((Real) i + 0.5) * dx[0];
640  Real y0 = ((Real) j + 0.5) * dx[1];
641 
642  // Final position at end of vector
643  x_pos_arr(i,j,k) = x0 + delta_x;
644  y_pos_arr(i,j,k) = y0 + delta_y;
645  Real z_bot_face = 0.25 * ( z_phys_arr(i ,j ,k) + z_phys_arr(i+1,j ,k)
646  + z_phys_arr(i ,j+1,k) + z_phys_arr(i+1,j+1,k) );
647  z_pos_arr(i,j,k) = z_bot_face + delta_z;
648 
649  // Destination position must be contained on the current process!
650  Real pos[] = {x_pos_arr(i,j,k),y_pos_arr(i,j,k),0.5*dx[2]};
651  amrex::ignore_unused(pos);
652  AMREX_ASSERT_WITH_MESSAGE( grb.contains(&pos[0]),
653  "Query point outside of proc domain!");
654  });
655  }
656 }

Referenced by make_MOSTAverage_at_level().

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

◆ set_plane_normalization()

void MOSTAverage::set_plane_normalization ( const int &  lev)

Function to compute normalization for plane average.

323 {
324  // Cells per plane and temp avg storage
325  m_ncell_plane.resize(m_maxlev);
326  m_plane_average.resize(m_maxlev);
327 
328  // True domain not used for normalization
329  Box domain = m_geom[lev].Domain();
330 
331  // NOTE: Level 0 spans the whole domain, but finer
332  // levels do not have such a restriction.
333  // For now, use the bounding box of the boxArray
334  // for normalization, consistent with avg routine.
335 
336  // Bounded box of CC data used for normalization
337  Box bnd_bx = (m_fields[lev][2]->boxArray()).minimalBox();
338 
339  // NOTE: Bounding box must lie on the periodic boundaries
340  // in order to trip the is_per flag
341 
342  // Num components, plane avg, cells per plane
343  Array<int,AMREX_SPACEDIM> is_per = {0,0,0};
344  for (int idim(0); idim < AMREX_SPACEDIM-1; ++idim) {
345  if ( m_geom[lev].isPeriodic(idim) &&
346  bnd_bx.bigEnd(idim)==domain.bigEnd(idim) &&
347  bnd_bx.smallEnd(idim)==domain.smallEnd(idim) ) { is_per[idim] = 1; }
348  }
349 
350  m_ncell_plane[lev].resize(m_navg);
351  m_plane_average[lev].resize(m_navg);
352  for (int iavg(0); iavg < m_navg; ++iavg) {
353  // Convert bnd_bx to current index type
354  IndexType ixt = m_averages[lev][iavg]->boxArray().ixType();
355  bnd_bx.convert(ixt);
356  IntVect bnd_bx_lo(bnd_bx.loVect());
357  IntVect bnd_bx_hi(bnd_bx.hiVect());
358 
359  m_plane_average[lev][iavg] = 0.0;
360 
361  m_ncell_plane[lev][iavg] = 1;
362  for (int idim(0); idim < AMREX_SPACEDIM; ++idim) {
363  if (idim != 2) {
364  if (ixt.nodeCentered(idim) && is_per[idim]) {
365  m_ncell_plane[lev][iavg] *= (bnd_bx_hi[idim] - bnd_bx_lo[idim]);
366  } else {
367  m_ncell_plane[lev][iavg] *= (bnd_bx_hi[idim] - bnd_bx_lo[idim] + 1);
368  }
369  }
370  } // idim
371  } // iavg
372 }

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
60  {m_ncell_region = (2 * m_radius + 1) * (2 * m_radius + 1) * (2 * m_radius + 1);}

Referenced by make_MOSTAverage_at_level().

Here is the caller graph for this function:

◆ set_rotated_fields()

void MOSTAverage::set_rotated_fields ( const int &  lev)

Function to set the rotated velocities.

264 {
265  // Peel back the level
266  auto& fields = m_fields[lev];
267  auto& rot_fields = m_rot_fields[lev];
268  auto z_phys_nd = m_z_phys_nd[lev];
269 
270  // Inverse grid size
271  const auto dxInv = m_geom[lev].InvCellSizeArray();
272 
273  // Single MFIter over CC data
274  int imf_cc = 2;
275 
276  // Populate rotated U & V for terrain
277 #ifdef _OPENMP
278 #pragma omp parallel if (Gpu::notInLaunchRegion())
279 #endif
280  for (MFIter mfi(*fields[imf_cc], TileNoZ()); mfi.isValid(); ++mfi) {
281  Box ubx = mfi.tilebox(IntVect(1,0,0));
282  Box vbx = mfi.tilebox(IntVect(0,1,0));
283 
284  const Array4<const Real>& z_phys_arr = z_phys_nd->const_array(mfi);
285 
286  const Array4<const Real>& u_arr = fields[0]->const_array(mfi);
287  const Array4<const Real>& v_arr = fields[1]->const_array(mfi);
288  const Array4<const Real>& w_arr = fields[4]->const_array(mfi);
289 
290  const Array4<Real>& u_rot_arr = rot_fields[0]->array(mfi);
291  const Array4<Real>& v_rot_arr = rot_fields[1]->array(mfi);
292 
293  // U rotated magnitude
294  ParallelFor(ubx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
295  {
296  // Elements of first tangent vector
297  Real met_h_xi = Compute_h_xi_AtIface(i,j,k,dxInv,z_phys_arr);
298  u_rot_arr(i,j,k) = (u_arr(i,j,k) + met_h_xi*w_arr(i,j,k))
299  / std::sqrt(met_h_xi*met_h_xi + 1.0);
300  });
301 
302  // V rotated magnitude
303  ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
304  {
305  // Elements of second tangent vector
306  Real met_h_eta = Compute_h_eta_AtJface(i,j,k,dxInv,z_phys_arr);
307  v_rot_arr(i,j,k) = (v_arr(i,j,k) + met_h_eta*w_arr(i,j,k))
308  / std::sqrt(met_h_eta*met_h_eta + 1.0);
309  });
310  }
311 
312  // Direct copy of other scalar variables
313  MultiFab::Copy(*rot_fields[2],*fields[2],0,0,1,rot_fields[2]->nGrowVect());
314  if (fields[3]) MultiFab::Copy(*rot_fields[3],*fields[3],0,0,1,rot_fields[3]->nGrowVect());
315 }
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:108
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:165

Referenced by compute_averages().

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

◆ set_z_positions_T()

void MOSTAverage::set_z_positions_T ( const int &  lev)

Function to set positions with terrain and e_z vector.

562 {
563  ParmParse pp(m_pp_prefix);
564  pp.get("most.zref",m_zref);
565 
566  // Capture for device
567  Real d_zref = m_zref;
568 
569  RealVect base;
570  const auto dx = m_geom[lev].CellSizeArray();
571  IntVect ng = m_x_pos[lev]->nGrowVect(); ng[2]=0;
572  for (MFIter mfi(*m_x_pos[lev], TileNoZ()); mfi.isValid(); ++mfi) {
573  Box npbx = mfi.tilebox(IntVect(1,1,0),IntVect(1,1,0));
574  Box gpbx = mfi.growntilebox(ng);
575  RealBox grb{gpbx,dx.data(),base.dataPtr()};
576 
577  const auto z_phys_arr = m_z_phys_nd[lev]->const_array(mfi);
578  auto x_pos_arr = m_x_pos[lev]->array(mfi);
579  auto y_pos_arr = m_y_pos[lev]->array(mfi);
580  auto z_pos_arr = m_z_pos[lev]->array(mfi);
581  ParallelFor(npbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
582  {
583  // Final position at end of vector
584  x_pos_arr(i,j,k) = ((Real) i + 0.5) * dx[0];
585  y_pos_arr(i,j,k) = ((Real) j + 0.5) * dx[1];
586  Real z_bot_face = 0.25 * ( z_phys_arr(i ,j ,k) + z_phys_arr(i+1,j ,k)
587  + z_phys_arr(i ,j+1,k) + z_phys_arr(i+1,j+1,k) );
588  z_pos_arr(i,j,k) = z_bot_face + d_zref;
589 
590  // Destination position must be contained on the current process!
591  Real pos[] = {x_pos_arr(i,j,k),y_pos_arr(i,j,k),0.5*dx[2]};
592  amrex::ignore_unused(pos);
593  AMREX_ASSERT_WITH_MESSAGE( grb.contains(&pos[0]),
594  "Query point outside of proc domain!");
595  });
596  }
597 }

Referenced by make_MOSTAverage_at_level().

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

◆ trilinear_interp_T()

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

Function to compute trilinear interpolation with terrain.

Parameters
[in]xpX-position
[in]ypY-position
[out]interp_valsValues interpolated
[in]interp_arrayArray to interpolate on
[in]z_arrPhysical heights
[in]ploProblem lower bounds
[in]dxiInverse cell size array
[in]interp_compNumber of components to interpolate
127  {
128  // Search to get z/k
129  bool found = false;
130  amrex::Real zval= 0.0;
131  int kmax = ubound(z_arr).z;
132  int i_new = (int) (xp * dxi[0] - 0.5);
133  int j_new = (int) (yp * dxi[1] - 0.5);
134  amrex::Real z_target = zp;
135  for (int lk(0); lk<kmax; ++lk) {
136  amrex::Real z_lo = 0.25 * ( z_arr(i_new,j_new ,lk ) + z_arr(i_new+1,j_new ,lk )
137  + z_arr(i_new,j_new+1,lk ) + z_arr(i_new+1,j_new+1,lk ) );
138  amrex::Real z_hi = 0.25 * ( z_arr(i_new,j_new ,lk+1) + z_arr(i_new+1,j_new ,lk+1)
139  + z_arr(i_new,j_new+1,lk+1) + z_arr(i_new+1,j_new+1,lk+1) );
140  if (z_target > z_lo && z_target < z_hi){
141  found = true;
142  zval = (amrex::Real) lk + ((z_target - z_lo) / (z_hi - z_lo)) + 0.5;
143  break;
144  }
145  }
146 
147  amrex::ignore_unused(found);
148  AMREX_ASSERT_WITH_MESSAGE(found, "MOSTAverage: Height above terrain not found, try increasing z_ref!");
149 
150  const amrex::RealVect lx((xp - plo[0])*dxi[0] + 0.5,
151  (yp - plo[1])*dxi[1] + 0.5,
152  zval);
153 
154  const amrex::IntVect ijk = lx.floor();
155 
156  int i = ijk[0]; int j = ijk[1]; int k = ijk[2];
157 
158  // Weights
159  const amrex::RealVect sx_hi = lx - ijk;
160  const amrex::RealVect sx_lo = 1.0 - sx_hi;
161 
162  for (int n = 0; n < interp_comp; n++)
163  interp_vals[n] = sx_lo[0]*sx_lo[1]*sx_lo[2]*interp_array(i-1, j-1, k-1,n) +
164  sx_lo[0]*sx_lo[1]*sx_hi[2]*interp_array(i-1, j-1, k ,n) +
165  sx_lo[0]*sx_hi[1]*sx_lo[2]*interp_array(i-1, j , k-1,n) +
166  sx_lo[0]*sx_hi[1]*sx_hi[2]*interp_array(i-1, j , k ,n) +
167  sx_hi[0]*sx_lo[1]*sx_lo[2]*interp_array(i , j-1, k-1,n) +
168  sx_hi[0]*sx_lo[1]*sx_hi[2]*interp_array(i , j-1, k ,n) +
169  sx_hi[0]*sx_hi[1]*sx_lo[2]*interp_array(i , j , k-1,n) +
170  sx_hi[0]*sx_hi[1]*sx_hi[2]*interp_array(i , j , k ,n);
171  }

Referenced by compute_plane_averages(), and compute_region_averages().

Here is the caller graph for this function:

◆ update_field_ptrs()

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

Function to reset the pointers to field variables.

Parameters
[in]levCurrent level
[in]vars_oldConserved variables at each level
[in]Theta_primPrimitive theta component at each level
249 {
250  m_fields[lev][0] = &vars_old[lev][Vars::xvel];
251  m_fields[lev][1] = &vars_old[lev][Vars::yvel];
252  m_fields[lev][2] = Theta_prim[lev].get();
253  m_fields[lev][3] = Qv_prim[lev].get();
254  m_fields[lev][4] = Qr_prim[lev].get();
255  m_fields[lev][5] = &vars_old[lev][Vars::zvel];
256 }

Referenced by ABLMost::update_mac_ptrs().

Here is the caller graph for this function:

◆ write_averages()

void MOSTAverage::write_averages ( const int &  lev)

Function to write averages to text file.

Parameters
[in]levCurrent level
1453 {
1454  // Peel back the level
1455  auto& averages = m_averages[lev];
1456 
1457  int navg = m_navg - 1;
1458 
1459  std::ofstream ofile;
1460  ofile.open ("MOST_averages.txt");
1461  ofile << "Averages computed via MOSTAverages class:\n";
1462 
1463  for (MFIter mfi(*averages[navg], TileNoZ()); mfi.isValid(); ++mfi) {
1464  Box bx = mfi.tilebox(); bx.setBig(2,0);
1465  int il = bx.smallEnd(0); int iu = bx.bigEnd(0);
1466  int jl = bx.smallEnd(1); int ju = bx.bigEnd(1);
1467 
1468  for (int j(jl); j <= ju; ++j) {
1469  for (int i(il); i <= iu; ++i) {
1470  ofile << "(I,J): " << "(" << i << "," << j << ")" << "\n";
1471  int k = 0;
1472  for (int iavg(0); iavg <= navg; ++iavg) {
1473  auto mf_arr = averages[iavg]->array(mfi);
1474  ofile << "iavg val: "
1475  << iavg << ' '
1476  << mf_arr(i,j,k) << "\n";
1477  }
1478  ofile << "\n";
1479  }
1480  }
1481  }
1482  ofile.close();
1483 }
Here is the call graph for this function:

◆ write_k_indices()

void MOSTAverage::write_k_indices ( const int &  lev)

Function to write the K indices to text file.

Parameters
[in]levCurrent level
1334 {
1335  // Peel back the level
1336  auto& averages = m_averages[lev];
1337  auto& k_indx = m_k_indx[lev];
1338 
1339  int navg = m_navg - 1;
1340 
1341  std::ofstream ofile;
1342  ofile.open ("MOST_k_indices.txt");
1343  ofile << "K indices used to compute averages via MOSTAverages class:\n";
1344 
1345  for (MFIter mfi(*averages[navg], TileNoZ()); mfi.isValid(); ++mfi) {
1346  Box bx = mfi.tilebox(); bx.setBig(2,0);
1347  int il = bx.smallEnd(0); int iu = bx.bigEnd(0);
1348  int jl = bx.smallEnd(1); int ju = bx.bigEnd(1);
1349 
1350  auto k_arr = k_indx->array(mfi);
1351 
1352  for (int j(jl); j <= ju; ++j) {
1353  for (int i(il); i <= iu; ++i) {
1354  ofile << "(I,J): " << "(" << i << "," << j << ")" << "\n";
1355  int k = 0;
1356  ofile << "K_ind: "
1357  << k_arr(i,j,k) << "\n";
1358  ofile << "\n";
1359  }
1360  }
1361  }
1362  ofile.close();
1363 }
Here is the call graph for this function:

◆ write_norm_indices()

void MOSTAverage::write_norm_indices ( const int &  lev)

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

Parameters
[in]levCurrent level
1373 {
1374  // Peel back the level
1375  auto& averages = m_averages[lev];
1376  auto& k_indx = m_k_indx[lev];
1377  auto& j_indx = m_j_indx[lev];
1378  auto& i_indx = m_i_indx[lev];
1379 
1380  int navg = m_navg - 1;
1381 
1382  std::ofstream ofile;
1383  ofile.open ("MOST_ijk_indices.txt");
1384  ofile << "IJK indices used to compute averages via MOSTAverages class:\n";
1385 
1386  for (MFIter mfi(*averages[navg], TileNoZ()); mfi.isValid(); ++mfi) {
1387  Box bx = mfi.tilebox(); bx.setBig(2,0);
1388  int il = bx.smallEnd(0); int iu = bx.bigEnd(0);
1389  int jl = bx.smallEnd(1); int ju = bx.bigEnd(1);
1390 
1391  auto k_arr = k_indx->array(mfi);
1392  auto j_arr = j_indx ? j_indx->array(mfi) : Array4<int> {};
1393  auto i_arr = i_indx ? i_indx->array(mfi) : Array4<int> {};
1394 
1395  for (int j(jl); j <= ju; ++j) {
1396  for (int i(il); i <= iu; ++i) {
1397  ofile << "(I1,J1,K1): " << "(" << i << "," << j << "," << 0 << ")" << "\n";
1398 
1399  int k = 0;
1400  int km = k_arr(i,j,k);
1401  int jm = j_arr ? j_arr(i,j,k) : j;
1402  int im = i_arr ? i_arr(i,j,k) : i;
1403 
1404  ofile << "(I2,J2,K2): "
1405  << "(" << im << "," << jm << "," << km << ")" << "\n";
1406  ofile << "\n";
1407  }
1408  }
1409  }
1410  ofile.close();
1411 }
Here is the call graph for this function:

◆ write_xz_positions()

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

Function to write X & Z positions to text file.

Parameters
[in]levCurrent level
[in]jIndex in y-dir
1423 {
1424  // Peel back the level
1425  auto& x_pos_mf = m_x_pos[lev];
1426  auto& z_pos_mf = m_z_pos[lev];
1427 
1428  std::ofstream ofile;
1429  ofile.open ("MOST_xz_positions.txt");
1430 
1431  for (MFIter mfi(*x_pos_mf, TileNoZ()); mfi.isValid(); ++mfi) {
1432  Box bx = mfi.tilebox(); bx.setBig(2,0);
1433  int il = bx.smallEnd(0); int iu = bx.bigEnd(0);
1434 
1435  auto x_pos_arr = x_pos_mf->array(mfi);
1436  auto z_pos_arr = z_pos_mf->array(mfi);
1437 
1438  int k = 0;
1439  for (int i(il); i <= iu; ++i)
1440  ofile << x_pos_arr(i,j,k) << ' ' << z_pos_arr(i,j,k) << "\n";
1441  }
1442  ofile.close();
1443 }
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_plane_averages(), compute_region_averages(), get_average(), make_MOSTAverage_at_level(), set_plane_normalization(), write_averages(), write_k_indices(), and write_norm_indices().

◆ m_fact_new

amrex::Real MOSTAverage::m_fact_new
protected

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

◆ m_fact_old

amrex::Real MOSTAverage::m_fact_old
protected

Referenced by 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_plane_averages(), compute_region_averages(), make_MOSTAverage_at_level(), set_plane_normalization(), set_rotated_fields(), and update_field_ptrs().

◆ m_geom

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

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

◆ m_i_indx

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

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

◆ m_interp

bool MOSTAverage::m_interp {false}
protected

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

◆ m_j_indx

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

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

◆ m_k_in

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

Referenced by set_k_indices_N(), and set_k_indices_T().

◆ m_k_indx

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

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

◆ m_maxlev

int MOSTAverage::m_maxlev {0}
protected

Referenced by make_MOSTAverage_at_level(), and set_plane_normalization().

◆ m_navg

int MOSTAverage::m_navg {6}
protected

Referenced by compute_plane_averages(), compute_region_averages(), make_MOSTAverage_at_level(), set_plane_normalization(), write_averages(), write_k_indices(), and write_norm_indices().

◆ m_ncell_plane

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

Referenced by compute_plane_averages(), and set_plane_normalization().

◆ m_ncell_region

int MOSTAverage::m_ncell_region {1}
protected

Referenced by compute_region_averages(), and set_region_normalization().

◆ m_norm_vec

bool MOSTAverage::m_norm_vec {false}
protected

Referenced by make_MOSTAverage_at_level().

◆ m_nvar

int MOSTAverage::m_nvar {6}
protected

Referenced by make_MOSTAverage_at_level().

◆ m_plane_average

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

Referenced by compute_plane_averages(), 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(), and set_z_positions_T().

◆ m_radius

int MOSTAverage::m_radius {0}
protected

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

◆ m_rot_fields

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

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

◆ m_rotate

bool MOSTAverage::m_rotate {false}
protected

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

◆ m_t_avg

bool MOSTAverage::m_t_avg {false}
protected

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

◆ m_t_init

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

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

◆ m_terrain_type

TerrainType MOSTAverage::m_terrain_type
protected

Referenced by make_MOSTAverage_at_level().

◆ m_time_window

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

Referenced by make_MOSTAverage_at_level().

◆ m_Vsg

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

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

◆ m_x_pos

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

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

◆ m_y_pos

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

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

◆ m_z_phys_nd

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

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

◆ m_z_pos

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

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

◆ m_zref

amrex::Real MOSTAverage::m_zref {10.0}
protected

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

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