#include <ERF_MOSTAverage.H>

Collaboration diagram for MOSTAverage:

[legend]

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

	~MOSTAverage ()

	MOSTAverage (MOSTAverage &&) noexcept=default

MOSTAverage &	operator= (MOSTAverage &&other) noexcept=delete

	MOSTAverage (const MOSTAverage &other)=delete

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

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

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

void	set_rotated_fields (const int &lev)

void	set_plane_normalization (const int &lev)

void	set_eb_normalization (const int &lev)

void	set_region_normalization (const int &)

void	set_k_indices_N (const int &lev)

void	set_k_indices_T (const int &lev)

void	set_norm_indices_T (const int &lev)

void	set_z_positions_T (const int &lev)

void	set_z_positions_EB (const int &lev)

void	set_norm_positions_T (const int &lev)

void	compute_averages (const int &lev)

void	compute_plane_averages (const int &lev)

void	compute_region_averages (const int &lev)

void	compute_eb_averages (const int &lev)

void	write_k_indices (const int &lev)

void	write_norm_indices (const int &lev)

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

void	write_averages (const int &lev)

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

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

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

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

Protected Attributes
const amrex::Vector< amrex::Geometry >	m_geom

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

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

std::string	m_pp_prefix

MeshType	m_mesh_type

TerrainType	m_terrain_type

int	m_nvar {6}

int	m_navg {6}

int	m_maxlev {0}

int	m_policy {0}

bool	m_rotate {false}

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

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

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

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

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

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

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

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

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

amrex::Vector< amrex::Vector< int > >	m_ncell_plane

amrex::Vector< amrex::Vector< amrex::Real > >	m_plane_average

int	m_radius {0}

int	m_ncell_region {1}

amrex::Vector< int >	m_k_in

bool	m_interp {false}

bool	m_norm_vec {false}

amrex::Vector< const eb_ * >	m_eb_vec

amrex::Vector< amrex::Vector< amrex::Real > >	m_total_bndry_area

bool	m_t_avg {false}

amrex::Vector< int >	m_t_init

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)

Constructor & Destructor Documentation

◆ MOSTAverage() [1/3]

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

explicit

Construct the MOST averaging helper.

Parameters

[in]	geom	geometry for all AMR levels
[in]	has_zphys	whether nodal physical-height data are available
[in]	a_pp_prefix	ParmParse prefix
[in]	m_mesh_type	mesh type
[in]	m_terrain_type	terrain representation
[in]	eb_vec	optional embedded-boundary geometry data

◆ ~MOSTAverage()

MOSTAverage::~MOSTAverage ( )

inline

Destroy the MOST averaging helper.

40 {}

◆ MOSTAverage() [2/3]

MOSTAverage::MOSTAverage ( MOSTAverage && )

defaultnoexcept

Default move constructor.

◆ MOSTAverage() [3/3]

MOSTAverage::MOSTAverage ( const MOSTAverage & other )

delete

Deleted copy constructor.

Parameters

[in] other source object

Member Function Documentation

◆ compute_averages()

void MOSTAverage::compute_averages ( const int & lev )

Driver for the selected average policy.

Parameters

[in] lev level index

Function to call the type of average computation.

Parameters

[in] lev Current level

 {
     if (m_rotate) set_rotated_fields(lev);
  
     switch(m_policy) {
     case 0: // Standard plane average
         compute_plane_averages(lev);
         break;
     case 1: // Local region/point
         compute_region_averages(lev);
         break;
     case 2: // EB Terrain average
         compute_eb_averages(lev);
         break;
     default:
         AMREX_ALWAYS_ASSERT_WITH_MESSAGE(false, "Unknown policy for MOSTAverage!");
     }
  
     // We have initialized the averages
     if (m_t_avg) m_t_init[lev] = 1;
 }

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] lev level index

Function to compute average over an EB surface.

Parameters

[in] lev Current level

 {
     AMREX_ALWAYS_ASSERT(m_eb_vec[lev] != nullptr);
  
     // Peel back the level
     auto& fields      = m_fields[lev];
     auto& averages    = m_averages[lev];
     auto& plane_average = m_plane_average[lev];
  
     // Get EB data
     const auto& cc_flags = m_eb_vec[lev]->get_const_factory()->getMultiEBCellFlagFab();
     auto cc_afrac = m_eb_vec[lev]->get_const_factory()->getAreaFrac();
     const auto& cc_bnorm = m_eb_vec[lev]->get_const_factory()->getBndryNormal();
     const auto& u_vfrac = m_eb_vec[lev]->get_u_const_factory()->getVolFrac();
     const auto& v_vfrac = m_eb_vec[lev]->get_v_const_factory()->getVolFrac();
     const auto& w_vfrac = m_eb_vec[lev]->get_w_const_factory()->getVolFrac();
  
     // Get geometry for cell sizes
     auto const& dx_arr = m_geom[lev].CellSizeArray();
     Real dx = dx_arr[0];
     Real dy = dx_arr[1];
     Real dz = dx_arr[2];
  
     // Set factors for time averaging
     Real d_fact_new, d_fact_old;
     if (m_t_avg && m_t_init[lev]) {
         d_fact_new = m_fact_new;
         d_fact_old = m_fact_old;
     } else {
         d_fact_new = one;
         d_fact_old = zero;
     }
  
     // GPU array to accumulate averages into
     Gpu::DeviceVector<Real> pavg(plane_average.size(), zero);
     Real* plane_avg = pavg.data();
  
     // Vectors for normalization and buffer storage
     Vector<Real> denom(plane_average.size(),zero);
     Vector<Real> val_old(plane_average.size(),zero);
  
     //
     //----------------------------------------------------------
     // Averages for U, V, and tangential velocity (in local coordinate)
     //----------------------------------------------------------
     //
     {
         denom[0]   = one / m_total_bndry_area[lev][0];
         val_old[0] = plane_average[0]*d_fact_old;
         denom[1]   = one / m_total_bndry_area[lev][1];
         val_old[1] = plane_average[1]*d_fact_old;
  
         int iavg = m_navg - 1;  // Tangential velocity magnitude
         denom[iavg]   = one / m_total_bndry_area[lev][iavg];
         val_old[iavg] = plane_average[iavg]*d_fact_old;
  
         const Real Vsg = m_Vsg[lev]; // Subgrid scale velocity
  
 #ifdef _OPENMP
 #pragma omp parallel if (Gpu::notInLaunchRegion())
 #endif
         for (MFIter mfi(*fields[2], TileNoZ()); mfi.isValid(); ++mfi) {
             const auto& flag = cc_flags[mfi];
  
             // Skip boxes that are not singlevalued (MultiCutFab only has data for singlevalued boxes)
             if (flag.getType() != FabType::singlevalued) continue;
  
             Box bx = mfi.tilebox();  // Full 3D box
  
             // Get EB arrays
             auto const flag_arr = flag.const_array();
             auto const afrac_x = cc_afrac[0]->const_array(mfi);
             auto const afrac_y = cc_afrac[1]->const_array(mfi);
             auto const afrac_z = cc_afrac[2]->const_array(mfi);
             auto const bnorm_arr = cc_bnorm.const_array(mfi);
             auto const u_vf_arr = u_vfrac.const_array(mfi);
             auto const v_vf_arr = v_vfrac.const_array(mfi);
             auto const w_vf_arr = w_vfrac.const_array(mfi);
  
             // Get velocity arrays
             auto const u_arr = fields[0]->const_array(mfi);
             auto const v_arr = fields[1]->const_array(mfi);
             auto const w_arr = fields[5]->const_array(mfi);
  
             ParallelFor(Gpu::KernelInfo().setReduction(true), bx, [=]
             AMREX_GPU_DEVICE(int i, int j, int k, Gpu::Handler const& handler) noexcept
             {
                 // Area-weighted averaging over cut cells at any k
                 if (flag_arr(i,j,k).isSingleValued()) {
                     // Compute area from face-centered area fractions
                     Real axm = afrac_x(i  ,j  ,k  );
                     Real axp = afrac_x(i+1,j  ,k  );
                     Real aym = afrac_y(i  ,j  ,k  );
                     Real ayp = afrac_y(i  ,j+1,k  );
                     Real azm = afrac_z(i  ,j  ,k  );
                     Real azp = afrac_z(i  ,j  ,k+1);
  
                     Real adx = (axm - axp) * dy * dz;
                     Real ady = (aym - ayp) * dx * dz;
                     Real adz = (azm - azp) * dx * dy;
  
                     Real area = std::sqrt(adx*adx + ady*ady + adz*adz);
  
                     // Volume-weighted interpolation of velocities to cell center
                     Real vf_u_lo = u_vf_arr(i,j,k);
                     Real vf_u_hi = u_vf_arr(i+1,j,k);
                     Real sum_vf_u = vf_u_lo + vf_u_hi;
                     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;
  
                     Real vf_v_lo = v_vf_arr(i,j,k);
                     Real vf_v_hi = v_vf_arr(i,j+1,k);
                     Real sum_vf_v = vf_v_lo + vf_v_hi;
                     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;
  
                     Real vf_w_lo = w_vf_arr(i,j,k);
                     Real vf_w_hi = w_vf_arr(i,j,k+1);
                     Real sum_vf_w = vf_w_lo + vf_w_hi;
                     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;
  
                     // Get normal vector components
                     Real nx = bnorm_arr(i,j,k,0);
                     Real ny = bnorm_arr(i,j,k,1);
                     Real nz = bnorm_arr(i,j,k,2);
  
                     // Compute tangential velocity components
                     Real v_dot_n = u_cc*nx + v_cc*ny + w_cc*nz;
                     Real u_tangent = u_cc - v_dot_n * nx;
                     Real v_tangent = v_cc - v_dot_n * ny;
                     Real mag = std::sqrt(u_tangent*u_tangent + v_tangent*v_tangent + Vsg*Vsg);
  
                     // Area-weighted sum
                     Real val_u = u_tangent * area;
                     Real val_v = v_tangent * area;
                     Real val_mag = mag * area;
  
                     Gpu::deviceReduceSum(&plane_avg[0], val_u, handler);
                     Gpu::deviceReduceSum(&plane_avg[1], val_v, handler);
                     Gpu::deviceReduceSum(&plane_avg[iavg], val_mag, handler);
                 }
             });
         }
     }
  
     //
     //----------------------------------------------------------
     // Averages for T,Qv (cell-centered scalars)
     //----------------------------------------------------------
     //
     for (int imf(2); imf < 4; ++imf) {
  
         // Continue if no valid Qv pointer
         if (!fields[imf]) continue;
  
         denom[imf]   = one / m_total_bndry_area[lev][imf];
         val_old[imf] = plane_average[imf]*d_fact_old;
  
 #ifdef _OPENMP
 #pragma omp parallel if (Gpu::notInLaunchRegion())
 #endif
         for (MFIter mfi(*fields[imf], TileNoZ()); mfi.isValid(); ++mfi) {
             const auto& flag = cc_flags[mfi];
  
             // Skip boxes that are not singlevalued (MultiCutFab only has data for singlevalued boxes)
             if (flag.getType() != FabType::singlevalued) continue;
  
             Box bx = mfi.tilebox();  // Full 3D box
             auto const flag_arr = flag.const_array();
             auto const afrac_x = cc_afrac[0]->const_array(mfi);
             auto const afrac_y = cc_afrac[1]->const_array(mfi);
             auto const afrac_z = cc_afrac[2]->const_array(mfi);
             auto const mf_arr = fields[imf]->const_array(mfi);
  
             ParallelFor(Gpu::KernelInfo().setReduction(true), bx, [=]
             AMREX_GPU_DEVICE(int i, int j, int k, Gpu::Handler const& handler) noexcept
             {
                 // Area-weighted averaging over cut cells at any k
                 if (flag_arr(i,j,k).isSingleValued()) {
                     // Compute area from face-centered area fractions
                     Real axm = afrac_x(i  ,j  ,k  );
                     Real axp = afrac_x(i+1,j  ,k  );
                     Real aym = afrac_y(i  ,j  ,k  );
                     Real ayp = afrac_y(i  ,j+1,k  );
                     Real azm = afrac_z(i  ,j  ,k  );
                     Real azp = afrac_z(i  ,j  ,k+1);
  
                     Real adx = (axm - axp) * dy * dz;
                     Real ady = (aym - ayp) * dx * dz;
                     Real adz = (azm - azp) * dx * dy;
  
                     Real area = std::sqrt(adx*adx + ady*ady + adz*adz);
  
                     Real val = mf_arr(i,j,k) * area;
                     Gpu::deviceReduceSum(&plane_avg[imf], val, handler);
                 }
             });
         }
     }
  
     //
     //------------------------------------------------------------------------
     // Averages for virtual potential temperature
     //------------------------------------------------------------------------
     //
     if (fields[3]) // We have water vapor
     {
         int iavg = 4;
         denom[iavg]   = one / m_total_bndry_area[lev][iavg];
         val_old[iavg] = plane_average[iavg]*d_fact_old;
  
 #ifdef _OPENMP
 #pragma omp parallel if (Gpu::notInLaunchRegion())
 #endif
         for (MFIter mfi(*fields[3], TileNoZ()); mfi.isValid(); ++mfi)
         {
             const auto& flag = cc_flags[mfi];
  
             // Skip boxes that are not singlevalued (MultiCutFab only has data for singlevalued boxes)
             if (flag.getType() != FabType::singlevalued) continue;
  
             Box bx = mfi.tilebox();  // Full 3D box
             auto const flag_arr = flag.const_array();
             auto const afrac_x = cc_afrac[0]->const_array(mfi);
             auto const afrac_y = cc_afrac[1]->const_array(mfi);
             auto const afrac_z = cc_afrac[2]->const_array(mfi);
  
             const Array4<Real const> T_mf_arr  = fields[2]->const_array(mfi);
             const Array4<Real const> qv_mf_arr = fields[3]->const_array(mfi);
             const Array4<Real const> qr_mf_arr = (fields[4]) ? fields[4]->const_array(mfi) :
                                                                 Array4<const Real> {};
  
             ParallelFor(Gpu::KernelInfo().setReduction(true), bx, [=]
             AMREX_GPU_DEVICE(int i, int j, int k, Gpu::Handler const& handler) noexcept
             {
                 if (flag_arr(i,j,k).isSingleValued()) {
                     // Compute area from face-centered area fractions
                     Real axm = afrac_x(i  ,j  ,k  );
                     Real axp = afrac_x(i+1,j  ,k  );
                     Real aym = afrac_y(i  ,j  ,k  );
                     Real ayp = afrac_y(i  ,j+1,k  );
                     Real azm = afrac_z(i  ,j  ,k  );
                     Real azp = afrac_z(i  ,j  ,k+1);
  
                     Real adx = (axm - axp) * dy * dz;
                     Real ady = (aym - ayp) * dx * dz;
                     Real adz = (azm - azp) * dx * dy;
  
                     Real area = std::sqrt(adx*adx + ady*ady + adz*adz);
  
                     Real vfac;
                     if (qr_mf_arr) {
                         // We also have liquid water
                         vfac = one + epsv*qv_mf_arr(i,j,k) - qr_mf_arr(i,j,k);
                     } else {
                         vfac = one + epsv*qv_mf_arr(i,j,k);
                     }
                     const Real val = T_mf_arr(i,j,k) * vfac * area;
                     Gpu::deviceReduceSum(&plane_avg[iavg], val, handler);
                 }
             });
         }
     }
     else // copy temperature
     {
         int iavg    = m_navg - 2;
         denom[iavg] = one / m_total_bndry_area[lev][iavg];
         // plane_avg[iavg] = plane_avg[2]
         Gpu::copy(Gpu::deviceToDevice, pavg.begin() + 2, pavg.begin() + 3,
                   pavg.begin() + iavg);
     }
  
     // Copy to host and sum across procs
     Gpu::copy(Gpu::deviceToHost, pavg.begin(), pavg.end(), plane_average.begin());
     ParallelDescriptor::ReduceRealSum(plane_average.data(), plane_average.size());
  
     // Normalize by total area and apply time averaging
     for (int iavg(0); iavg < m_navg; ++iavg){
         plane_average[iavg] *= denom[iavg]*d_fact_new;
         plane_average[iavg] += val_old[iavg];
         averages[iavg]->setVal(plane_average[iavg]);
     }
 }

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] lev level index

Function to compute average over a plane.

Parameters

[in] lev Current level

 {
     // Peel back the level
     auto& fields      = m_fields[lev];
     auto& rot_fields  = m_rot_fields[lev];
     auto& averages    = m_averages[lev];
     const auto & geom = m_geom[lev];
  
     auto& z_phys   = m_z_phys_nd[lev];
     auto& x_pos    = m_x_pos[lev];
     auto& y_pos    = m_y_pos[lev];
     auto& z_pos    = m_z_pos[lev];
  
     auto& i_indx   = m_i_indx[lev];
     auto& j_indx   = m_j_indx[lev];
     auto& k_indx   = m_k_indx[lev];
  
     auto& ncell_plane   = m_ncell_plane[lev];
     auto& plane_average = m_plane_average[lev];
  
     // Set factors for time averaging
     Real d_fact_new, d_fact_old;
     if (m_t_avg && m_t_init[lev]) {
         d_fact_new = m_fact_new;
         d_fact_old = m_fact_old;
     } else {
         d_fact_new = one;
         d_fact_old = zero;
     }
  
     int klo = m_geom[lev].Domain().smallEnd(2);
  
     // GPU array to accumulate averages into
     Gpu::DeviceVector<Real> pavg(plane_average.size(), zero);
     Real* plane_avg = pavg.data();
  
     // Vectors for normalization and buffer storage
     Vector<Real> denom(plane_average.size(),zero);
     Vector<Real> val_old(plane_average.size(),zero);
  
     //
     //----------------------------------------------------------
     // Averages over all the fields
     //----------------------------------------------------------
     //
     Box domain = geom.Domain();
  
     Array<int,AMREX_SPACEDIM> is_per = {0,0,0};
     for (int idim(0); idim < AMREX_SPACEDIM-1; ++idim) {
         if (geom.isPeriodic(idim)) is_per[idim] = 1;
     }
  
     // Averages for U,V,T,Qv (not Qc or W)
     for (int imf(0); imf < 4; ++imf) {
  
         // Continue if no valid Qv pointer
         if (!fields[imf]) continue;
  
         denom[imf]   = one / (Real)ncell_plane[imf];
         val_old[imf] = plane_average[imf]*d_fact_old;
  
 #ifdef _OPENMP
 #pragma omp parallel if (Gpu::notInLaunchRegion())
 #endif
         for (MFIter mfi(*fields[imf], TileNoZ()); mfi.isValid(); ++mfi) {
             Box vbx = mfi.validbox(); // This is the grid (not tile)
             Box pbx = mfi.tilebox();  // This is the tile (not grid)
  
             if (pbx.smallEnd(2) != klo) { continue; }
  
             // Make planar since mfiter is over fields
             pbx.makeSlab(2,klo);
  
             // Avoid double counting nodal data by changing the high end when we are
             //     at the high side of the grid (not just of the tile)
             IndexType ixt = averages[imf]->boxArray().ixType();
             for (int idim(0); idim < AMREX_SPACEDIM-1; ++idim) {
                 if ( ixt.nodeCentered(idim)  && (pbx.bigEnd(idim) == vbx.bigEnd(idim)) ) {
                     int dom_hi = domain.bigEnd(idim)+1;
                     if (pbx.bigEnd(idim) < dom_hi || is_per[idim]) {
                         pbx.growHi(idim,-1);
                     }
                 }
             }
  
             auto mf_arr = (m_rotate) ? rot_fields[imf]->const_array(mfi) :
                                            fields[imf]->const_array(mfi);
  
             if (m_interp) {
                 const auto plo   = geom.ProbLoArray();
                 const auto dxInv = geom.InvCellSizeArray();
                 const auto z_phys_arr = z_phys->const_array(mfi);
                 auto x_pos_arr = x_pos->array(mfi);
                 auto y_pos_arr = y_pos->array(mfi);
                 auto z_pos_arr = z_pos->array(mfi);
                 ParallelFor(Gpu::KernelInfo().setReduction(true), pbx, [=]
                 AMREX_GPU_DEVICE(int i, int j, int , Gpu::Handler const& handler) noexcept
                 {
                     Real interp{0};
                     trilinear_interp_T(x_pos_arr(i,j,0), y_pos_arr(i,j,0), z_pos_arr(i,j,0),
                                        &interp, mf_arr, z_phys_arr, plo, dxInv, 1);
                     Real val = interp;
                     Gpu::deviceReduceSum(&plane_avg[imf], val, handler);
                 });
             } else {
                 auto k_arr  = k_indx->const_array(mfi);
                 auto j_arr  = j_indx ? j_indx->const_array(mfi) : Array4<const int> {};
                 auto i_arr  = i_indx ? i_indx->const_array(mfi) : Array4<const int> {};
                 ParallelFor(Gpu::KernelInfo().setReduction(true), pbx, [=]
                 AMREX_GPU_DEVICE(int i, int j, int , Gpu::Handler const& handler) noexcept
                 {
                     int mk = k_arr(i,j,0);
                     int mj = j_arr ? j_arr(i,j,0) : j;
                     int mi = i_arr ? i_arr(i,j,0) : i;
                     Real val = mf_arr(mi,mj,mk);
                     Gpu::deviceReduceSum(&plane_avg[imf], val, handler);
                 });
             }
         }
     }
  
     //
     //------------------------------------------------------------------------
     // Averages for virtual potential temperature
     // (This is cell-centered so we don't need to worry about double-counting)
     //------------------------------------------------------------------------
     //
     if (fields[3]) // We have water vapor
     {
         int iavg = 4;
         denom[iavg]   = one / (Real)ncell_plane[iavg];
         val_old[iavg] = plane_average[iavg]*d_fact_old;
  
 #ifdef _OPENMP
 #pragma omp parallel if (Gpu::notInLaunchRegion())
 #endif
         for (MFIter mfi(*fields[3], TileNoZ()); mfi.isValid(); ++mfi)
         {
             Box pbx = mfi.tilebox();
  
             if (pbx.smallEnd(2) != klo) { continue; }
  
             pbx.makeSlab(2,klo);
  
             const Array4<Real const>& T_mf_arr  = fields[2]->const_array(mfi);
             const Array4<Real const>& qv_mf_arr = fields[3]->const_array(mfi);
             const Array4<Real const>& qr_mf_arr = (fields[4]) ? fields[4]->const_array(mfi) :
                                                                 Array4<const Real> {};
  
             if (m_interp) {
                 const auto plo   = m_geom[lev].ProbLoArray();
                 const auto dxInv = m_geom[lev].InvCellSizeArray();
                 const auto z_phys_arr = z_phys->const_array(mfi);
                 auto x_pos_arr = x_pos->array(mfi);
                 auto y_pos_arr = y_pos->array(mfi);
                 auto z_pos_arr = z_pos->array(mfi);
                 ParallelFor(Gpu::KernelInfo().setReduction(true), pbx, [=]
                 AMREX_GPU_DEVICE(int i, int j, int , Gpu::Handler const& handler) noexcept
                 {
                     Real T_interp{0};
                     Real qv_interp{0};
                     trilinear_interp_T(x_pos_arr(i,j,0), y_pos_arr(i,j,0), z_pos_arr(i,j,0),
                                        &T_interp, T_mf_arr, z_phys_arr, plo, dxInv, 1);
                     trilinear_interp_T(x_pos_arr(i,j,0), y_pos_arr(i,j,0), z_pos_arr(i,j,0),
                                        &qv_interp, qv_mf_arr, z_phys_arr, plo, dxInv, 1);
                     Real vfac;
                     if (qr_mf_arr) {
                         // We also have liquid water
                         Real qr_interp{0};
                         trilinear_interp_T(x_pos_arr(i,j,0), y_pos_arr(i,j,0), z_pos_arr(i,j,0),
                                            &qr_interp, qr_mf_arr, z_phys_arr, plo, dxInv, 1);
                         vfac = one + epsv*qv_interp - qr_interp;
                     } else {
                         vfac = one + epsv*qv_interp;
                     }
                     const Real val = T_interp * vfac;
                     Gpu::deviceReduceSum(&plane_avg[iavg], val, handler);
                 });
             } else {
                 auto k_arr = k_indx->const_array(mfi);
                 auto j_arr = j_indx ? j_indx->const_array(mfi) : Array4<const int> {};
                 auto i_arr = i_indx ? i_indx->const_array(mfi) : Array4<const int> {};
                 ParallelFor(Gpu::KernelInfo().setReduction(true), pbx, [=]
                 AMREX_GPU_DEVICE(int i, int j, int , Gpu::Handler const& handler) noexcept
                 {
                     int mk = k_arr(i,j,0);
                     int mj = j_arr ? j_arr(i,j,0) : j;
                     int mi = i_arr ? i_arr(i,j,0) : i;
                     Real vfac;
                     if (qr_mf_arr) {
                         // We also have liquid water
                         vfac = one + epsv*qv_mf_arr(mi,mj,mk) - qr_mf_arr(mi,mj,mk);
                     } else {
                         vfac = one + epsv*qv_mf_arr(mi,mj,mk);
                     }
                     const Real val = T_mf_arr(mi,mj,mk) * vfac;
                     Gpu::deviceReduceSum(&plane_avg[iavg], val, handler);
                 });
             }
         }
     }
     else // copy temperature
     {
         int iavg    = m_navg - 2;
         denom[iavg] = one / (Real)ncell_plane[iavg];
         // plane_avg[iavg] = plane_avg[2]
         Gpu::copy(Gpu::deviceToDevice, pavg.begin() + 2, pavg.begin() + 3,
                   pavg.begin() + iavg);
     }
  
     //
     //------------------------------------------------------------------------
     // Averages for the tangential velocity magnitude
     // (This is cell-centered so we don't need to worry about double-counting)
     //------------------------------------------------------------------------
     //
     {
         int imf_cc = 2;
         int imf  = 0;
         int iavg = m_navg - 1;
         denom[iavg]   = one / (Real)ncell_plane[iavg];
         val_old[iavg] = plane_average[iavg]*d_fact_old;
  
         const Real Vsg = m_Vsg[lev];
  
 #ifdef _OPENMP
 #pragma omp parallel if (Gpu::notInLaunchRegion())
 #endif
         for (MFIter mfi(*fields[imf_cc], TileNoZ()); mfi.isValid(); ++mfi)
         {
             Box pbx = mfi.tilebox();
  
             if (pbx.smallEnd(2) != klo) { continue; }
  
             pbx.makeSlab(2,klo);
  
             // Last element is Umag and always cell centered
             auto u_mf_arr = (m_rotate) ? rot_fields[imf  ]->const_array(mfi) :
                                              fields[imf  ]->const_array(mfi);
             auto v_mf_arr = (m_rotate) ? rot_fields[imf+1]->const_array(mfi) :
                                              fields[imf+1]->const_array(mfi);
  
             if (m_interp) {
                 const auto plo   = m_geom[lev].ProbLoArray();
                 const auto dxInv = m_geom[lev].InvCellSizeArray();
                 const auto z_phys_arr = z_phys->const_array(mfi);
                 auto x_pos_arr = x_pos->array(mfi);
                 auto y_pos_arr = y_pos->array(mfi);
                 auto z_pos_arr = z_pos->array(mfi);
                 ParallelFor(Gpu::KernelInfo().setReduction(true), pbx, [=]
                 AMREX_GPU_DEVICE(int i, int j, int , Gpu::Handler const& handler) noexcept
                 {
                     Real u_interp{0};
                     Real v_interp{0};
                     trilinear_interp_T(x_pos_arr(i,j,0), y_pos_arr(i,j,0), z_pos_arr(i,j,0),
                                        &u_interp, u_mf_arr, z_phys_arr, plo, dxInv, 1);
                     trilinear_interp_T(x_pos_arr(i,j,0), y_pos_arr(i,j,0), z_pos_arr(i,j,0),
                                        &v_interp, v_mf_arr, z_phys_arr, plo, dxInv, 1);
                     const Real val = std::sqrt(u_interp*u_interp + v_interp*v_interp + Vsg*Vsg);
                     Gpu::deviceReduceSum(&plane_avg[iavg], val, handler);
                 });
             } else {
                 auto k_arr = k_indx->const_array(mfi);
                 auto j_arr = j_indx ? j_indx->const_array(mfi) : Array4<const int> {};
                 auto i_arr = i_indx ? i_indx->const_array(mfi) : Array4<const int> {};
                 ParallelFor(Gpu::KernelInfo().setReduction(true), pbx, [=]
                 AMREX_GPU_DEVICE(int i, int j, int , Gpu::Handler const& handler) noexcept
                 {
                     int mk = k_arr(i,j,0);
                     int mj = j_arr ? j_arr(i,j,0) : j;
                     int mi = i_arr ? i_arr(i,j,0) : i;
                     const Real u_val = myhalf * (u_mf_arr(mi,mj,mk) + u_mf_arr(mi+1,mj  ,mk));
                     const Real v_val = myhalf * (v_mf_arr(mi,mj,mk) + v_mf_arr(mi  ,mj+1,mk));
                     const Real val = std::sqrt(u_val*u_val + v_val*v_val + Vsg*Vsg);
                     Gpu::deviceReduceSum(&plane_avg[iavg], val, handler);
                 });
             }
         }
     }
  
     // Copy to host and sum across procs
     Gpu::copy(Gpu::deviceToHost, pavg.begin(), pavg.end(), plane_average.begin());
     ParallelDescriptor::ReduceRealSum(plane_average.data(), static_cast<int>(plane_average.size()));
  
     // No spatial variation with plane averages
     for (int iavg(0); iavg < m_navg; ++iavg){
         plane_average[iavg] *= denom[iavg]*d_fact_new;
         plane_average[iavg] += val_old[iavg];
         averages[iavg]->setVal(plane_average[iavg]);
     }
 }

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] lev level index

Function to compute average over local region.

Parameters

[in] lev Current level

 {
     // Peel back the level
     auto& fields      = m_fields[lev];
     auto& rot_fields  = m_rot_fields[lev];
     auto& averages    = m_averages[lev];
     const auto & geom = m_geom[lev];
  
     auto& z_phys   = m_z_phys_nd[lev];
     auto& x_pos    = m_x_pos[lev];
     auto& y_pos    = m_y_pos[lev];
     auto& z_pos    = m_z_pos[lev];
  
     auto& i_indx   = m_i_indx[lev];
     auto& j_indx   = m_j_indx[lev];
     auto& k_indx   = m_k_indx[lev];
  
     int klo = m_geom[lev].Domain().smallEnd(2);
  
     // Set factors for time averaging
     Real d_fact_new, d_fact_old;
     if (m_t_avg && m_t_init[lev]) {
         d_fact_new = m_fact_new;
         d_fact_old = m_fact_old;
     } else {
         d_fact_new = one;
         d_fact_old = zero;
     }
  
     // Number of cells contained in the local average
     const Real denom = one / (Real) m_ncell_region;
  
     // Capture radius for device
     int d_radius = m_radius;
  
     //
     //----------------------------------------------------------
     // Averages for U,V,T,Qv
     //----------------------------------------------------------
     //
     for (int imf(0); imf < 4; ++imf) {
  
         // Continue if no valid Qv pointer
         if (!fields[imf]) continue;
  
 #ifdef _OPENMP
 #pragma omp parallel if (Gpu::notInLaunchRegion())
 #endif
         for (MFIter mfi(*fields[imf], TileNoZ()); mfi.isValid(); ++mfi) {
             Box pbx = mfi.tilebox();
  
             if (pbx.smallEnd(2) != klo) { continue; }
  
             // Make planar since mfiter is over fields
             pbx.makeSlab(2,klo);
  
             auto mf_arr = (m_rotate) ? rot_fields[imf]->const_array(mfi) :
                                            fields[imf]->const_array(mfi);
             auto ma_arr = averages[imf]->array(mfi);
  
             if (m_interp) {
                 const auto plo   = geom.ProbLoArray();
                 const auto dx    = geom.CellSizeArray();
                 const auto dxInv = geom.InvCellSizeArray();
                 const auto z_phys_arr = z_phys->const_array(mfi);
                 auto x_pos_arr = x_pos->array(mfi);
                 auto y_pos_arr = y_pos->array(mfi);
                 auto z_pos_arr = z_pos->array(mfi);
                 ParallelFor(pbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
                 {
                     ma_arr(i,j,0) *= d_fact_old;
  
                     Real met_h_zeta = Compute_h_zeta_AtCellCenter(i,j,k,dxInv,z_phys_arr);
                     for (int lk(-d_radius); lk <= (d_radius); ++lk) {
                       for (int lj(-d_radius); lj <= (d_radius); ++lj) {
                         for (int li(-d_radius); li <= (d_radius); ++li) {
                             Real interp{0};
                             Real xp = x_pos_arr(i+li,j+lj,0);
                             Real yp = y_pos_arr(i+li,j+lj,0);
                             Real zp = z_pos_arr(i+li,j+lj,0) + met_h_zeta*lk*dx[2];
                             trilinear_interp_T(xp, yp, zp, &interp, mf_arr, z_phys_arr, plo, dxInv, 1);
                             Real val = denom * interp * d_fact_new;
                             ma_arr(i,j,0) += val;
                         }
                       }
                     }
                 });
             } else {
                 auto k_arr = k_indx->const_array(mfi);
                 auto j_arr = j_indx ? j_indx->const_array(mfi) : Array4<const int> {};
                 auto i_arr = i_indx ? i_indx->const_array(mfi) : Array4<const int> {};
                 ParallelFor(pbx, [=] AMREX_GPU_DEVICE(int i, int j, int ) noexcept
                 {
                     ma_arr(i,j,0) *= d_fact_old;
  
                     int mk = k_arr(i,j,0);
                     int mj = j_arr ? j_arr(i,j,0) : j;
                     int mi = i_arr ? i_arr(i,j,0) : i;
                     for (int lk(mk-d_radius); lk <= (mk+d_radius); ++lk) {
                       for (int lj(mj-d_radius); lj <= (mj+d_radius); ++lj) {
                         for (int li(mi-d_radius); li <= (mi+d_radius); ++li) {
                             Real val = denom * mf_arr(li, lj, lk) * d_fact_new;
                             ma_arr(i,j,0) += val;
                         }
                       }
                     }
                 });
             }
         } // MFiter
  
         // Fill interior ghost cells and any ghost cells outside a periodic domain
         //***********************************************************************************
         averages[imf]->FillBoundary(geom.periodicity());
  
     } // imf
  
     //
     //----------------------------------------------------------
     // Averages for virtual potential temperature
     //----------------------------------------------------------
     //
     if (fields[3]) // We have water vapor
     {
         int iavg = 4;
  
 #ifdef _OPENMP
 #pragma omp parallel if (Gpu::notInLaunchRegion())
 #endif
         for (MFIter mfi(*fields[3], TileNoZ()); mfi.isValid(); ++mfi) {
             Box pbx = mfi.tilebox();
  
             if (pbx.smallEnd(2) != klo) { continue; }
  
             pbx.makeSlab(2,klo);
  
             const Array4<Real const>& T_mf_arr  = fields[2]->const_array(mfi);
             const Array4<Real const>& qv_mf_arr = fields[3]->const_array(mfi);
             const Array4<Real const>& qr_mf_arr = (fields[4]) ? fields[4]->const_array(mfi) :
                                                                 Array4<const Real> {};
             auto ma_arr = averages[iavg]->array(mfi);
  
             if (m_interp) {
                 const auto plo   = geom.ProbLoArray();
                 const auto dx    = geom.CellSizeArray();
                 const auto dxInv = geom.InvCellSizeArray();
                 const auto z_phys_arr = z_phys->const_array(mfi);
                 auto x_pos_arr = x_pos->array(mfi);
                 auto y_pos_arr = y_pos->array(mfi);
                 auto z_pos_arr = z_pos->array(mfi);
                 ParallelFor(pbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
                 {
                     ma_arr(i,j,0) *= d_fact_old;
  
                     Real met_h_zeta = Compute_h_zeta_AtCellCenter(i,j,k,dxInv,z_phys_arr);
                     for (int lk(-d_radius); lk <= (d_radius); ++lk) {
                       for (int lj(-d_radius); lj <= (d_radius); ++lj) {
                         for (int li(-d_radius); li <= (d_radius); ++li) {
                             Real T_interp{0};
                             Real qv_interp{0};
                             Real xp = x_pos_arr(i+li,j+lj,0);
                             Real yp = y_pos_arr(i+li,j+lj,0);
                             Real zp = z_pos_arr(i+li,j+lj,0) + met_h_zeta*lk*dx[2];
                             trilinear_interp_T(xp, yp, zp, &T_interp,  T_mf_arr,  z_phys_arr, plo, dxInv, 1);
                             trilinear_interp_T(xp, yp, zp, &qv_interp, qv_mf_arr, z_phys_arr, plo, dxInv, 1);
                             Real vfac;
                             if (qr_mf_arr) {
                                 // We also have liquid water
                                 Real qr_interp{0};
                                 trilinear_interp_T(x_pos_arr(i,j,0), y_pos_arr(i,j,0), z_pos_arr(i,j,0),
                                                    &qr_interp, qr_mf_arr, z_phys_arr, plo, dxInv, 1);
                                 vfac = one + epsv*qv_interp - qr_interp;
                             } else {
                                 vfac = one + epsv*qv_interp;
                             }
                             const Real mag = T_interp * vfac;
                             const Real val = denom * mag * d_fact_new;
                             ma_arr(i,j,0) += val;
                         }
                       }
                     }
                 });
             } else {
                 auto k_arr = k_indx->const_array(mfi);
                 auto j_arr = j_indx ? j_indx->const_array(mfi) : Array4<const int> {};
                 auto i_arr = i_indx ? i_indx->const_array(mfi) : Array4<const int> {};
                 ParallelFor(pbx, [=] AMREX_GPU_DEVICE(int i, int j, int ) noexcept
                 {
                     ma_arr(i,j,0) *= d_fact_old;
  
                     int mk = k_arr(i,j,0);
                     int mj = j_arr ? j_arr(i,j,0) : j;
                     int mi = i_arr ? i_arr(i,j,0) : i;
                     for (int lk(mk-d_radius); lk <= (mk+d_radius); ++lk) {
                       for (int lj(mj-d_radius); lj <= (mj+d_radius); ++lj) {
                         for (int li(mi-d_radius); li <= (mi+d_radius); ++li) {
                             Real vfac;
                             if (qr_mf_arr) {
                                 // We also have liquid water
                                 vfac = one + epsv*qv_mf_arr(li,lj,lk) - qr_mf_arr(li,lj,lk);
                             } else {
                                 vfac = one + epsv*qv_mf_arr(li,lj,lk);
                             }
                             const Real mag = T_mf_arr(li,lj,lk) * vfac;
                             const Real val = denom * mag * d_fact_new;
                             ma_arr(i,j,0) += val;
                         }
                       }
                     }
                 });
             }
         } // MFiter
  
         // Fill interior ghost cells and any ghost cells outside a periodic domain
         //***********************************************************************************
         averages[iavg]->FillBoundary(geom.periodicity());
  
     }
     else // copy temperature
     {
         int iavg   = m_navg - 2;
         IntVect ng = averages[iavg]->nGrowVect();
         MultiFab::Copy(*(averages[iavg]),*(averages[2]),0,0,1,ng);
     }
  
     //
     //----------------------------------------------------------
     // Averages for the tangential velocity magnitude
     //----------------------------------------------------------
     //
     {
         int imf_cc = 2;
         int imf  = 0;
         int iavg = m_navg - 1;
  
         const Real Vsg = m_Vsg[lev];
  
 #ifdef _OPENMP
 #pragma omp parallel if (Gpu::notInLaunchRegion())
 #endif
         for (MFIter mfi(*fields[imf_cc], TileNoZ()); mfi.isValid(); ++mfi) {
             Box pbx = mfi.tilebox();
  
             if (pbx.smallEnd(2) != klo) { continue; }
  
             pbx.makeSlab(2,klo);
  
             auto u_mf_arr = (m_rotate) ? rot_fields[imf  ]->const_array(mfi) :
                                              fields[imf  ]->const_array(mfi);
             auto v_mf_arr = (m_rotate) ? rot_fields[imf+1]->const_array(mfi) :
                                              fields[imf+1]->const_array(mfi);
             auto ma_arr   = averages[iavg]->array(mfi);
  
             if (m_interp) {
                 const auto plo   = geom.ProbLoArray();
                 const auto dx    = geom.CellSizeArray();
                 const auto dxInv = geom.InvCellSizeArray();
                 const auto z_phys_arr = z_phys->const_array(mfi);
                 auto x_pos_arr = x_pos->array(mfi);
                 auto y_pos_arr = y_pos->array(mfi);
                 auto z_pos_arr = z_pos->array(mfi);
                 ParallelFor(pbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
                 {
                     ma_arr(i,j,0) *= d_fact_old;
  
                     Real met_h_zeta = Compute_h_zeta_AtCellCenter(i,j,k,dxInv,z_phys_arr);
                     for (int lk(-d_radius); lk <= (d_radius); ++lk) {
                       for (int lj(-d_radius); lj <= (d_radius); ++lj) {
                         for (int li(-d_radius); li <= (d_radius); ++li) {
                             Real u_interp{0};
                             Real v_interp{0};
                             Real xp = x_pos_arr(i+li,j+lj,0);
                             Real yp = y_pos_arr(i+li,j+lj,0);
                             Real zp = z_pos_arr(i+li,j+lj,0) + met_h_zeta*lk*dx[2];
                             trilinear_interp_T(xp, yp, zp, &u_interp, u_mf_arr, z_phys_arr, plo, dxInv, 1);
                             trilinear_interp_T(xp, yp, zp, &v_interp, v_mf_arr, z_phys_arr, plo, dxInv, 1);
                             const Real mag = std::sqrt(u_interp*u_interp + v_interp*v_interp + Vsg*Vsg);
                             Real val = denom * mag * d_fact_new;
                             ma_arr(i,j,0) += val;
                         }
                       }
                     }
                 });
             } else {
                 auto k_arr = k_indx->const_array(mfi);
                 auto j_arr = j_indx ? j_indx->const_array(mfi) : Array4<const int> {};
                 auto i_arr = i_indx ? i_indx->const_array(mfi) : Array4<const int> {};
                 ParallelFor(pbx, [=] AMREX_GPU_DEVICE(int i, int j, int ) noexcept
                 {
                     ma_arr(i,j,0) *= d_fact_old;
  
                     int mk = k_arr(i,j,0);
                     int mj = j_arr ? j_arr(i,j,0) : j;
                     int mi = i_arr ? i_arr(i,j,0) : i;
                     for (int lk(mk-d_radius); lk <= (mk+d_radius); ++lk) {
                       for (int lj(mj-d_radius); lj <= (mj+d_radius); ++lj) {
                         for (int li(mi-d_radius); li <= (mi+d_radius); ++li) {
                             const Real u_val = myhalf * (u_mf_arr(li,lj,lk) + u_mf_arr(li+1,lj  ,lk));
                             const Real v_val = myhalf * (v_mf_arr(li,lj,lk) + v_mf_arr(li  ,lj+1,lk));
                             const Real mag   = std::sqrt(u_val*u_val + v_val*v_val + Vsg*Vsg);
                             Real val = denom * mag * d_fact_new;
                             ma_arr(i,j,0) += val;
                         }
                       }
                     }
                 });
             }
         } // MFiter
  
         // Fill interior ghost cells and any ghost cells outside a periodic domain
         //***********************************************************************************
         averages[iavg]->FillBoundary(geom.periodicity());
  
     }
  
     // NOTE: Checking periodicity with the geom structure is not
     //       sufficient at higher levels. The BA may be contained
     //       within the domain and it's exterior ghost cells filled
     //       from interpolation; yet the domain BCs are periodic.
  
     // Need to fill ghost cells outside the domain if not periodic
     bool not_per_x = !(geom.periodicity().isPeriodic(0));
     bool not_per_y = !(geom.periodicity().isPeriodic(1));
     Box cc_bnd_bx  = (m_fields[lev][2]->boxArray()).minimalBox();
     Box domain     = geom.Domain();
     if (domain.contains(cc_bnd_bx) || (not_per_x || not_per_y)) {
         for (int iavg(0); iavg < m_navg; ++iavg) {
             IntVect ng = averages[iavg]->nGrowVect(); ng[2]=0;
  
             // NOTE:  Level 0 spans the whole domain, but finer
             //        levels do not have such a restriction.
             //        For now, use the bounding box of the boxArray.
  
             // NOTE2: The fields and averages have different indexing.
             //        The averages are: U/V/T/Qv/Tv/Umag
             //        The fields   are: U/V/T/Qv/Qr/W
             //        We clip iavg at 2 since all the remaining data is CC
  
             // Bounded box of CC data used for normalization
             int imf = min(iavg,2);
             Box bnd_bx = (fields[imf]->boxArray()).minimalBox();
 #ifdef _OPENMP
 #pragma omp parallel if (Gpu::notInLaunchRegion())
 #endif
             for (MFIter mfi(*fields[imf], TileNoZ()); mfi.isValid(); ++mfi) {
                 Box gpbx = mfi.growntilebox(ng);
  
                 if (gpbx.smallEnd(2) != klo) { continue; }
  
                 gpbx.makeSlab(2,klo);
  
                 if (bnd_bx.contains(gpbx)) continue;
  
                 auto ma_arr = averages[iavg]->array(mfi);
  
                 int i_lo = bnd_bx.smallEnd(0); int i_hi = bnd_bx.bigEnd(0);
                 int j_lo = bnd_bx.smallEnd(1); int j_hi = bnd_bx.bigEnd(1);
                 ParallelFor(gpbx, [=] AMREX_GPU_DEVICE(int i, int j, int ) noexcept
                 {
                     int li, lj;
                     li = i  < i_lo ? i_lo : i;
                     li = li > i_hi ? i_hi : li;
                     lj = j  < j_lo ? j_lo : j;
                     lj = lj > j_hi ? j_hi : lj;
  
                     ma_arr(i,j,0) = ma_arr(li,lj,0);
                 });
             } // MFiter
         } // iavg
     } // Not periodic
 }

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

Return one 2D average MultiFab.

Parameters

[in]	lev	level index
[in]	comp	average component index

235 { return m_averages[lev][comp].get(); }

Referenced by SurfaceLayer::get_mac_avg().

Here is the caller graph for this function:

◆ get_k_indices()

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

inline

Return the k-index iMultiFab.

Parameters

[in] lev level index

249 { return m_k_indx[lev].get(); }

◆ 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] lev level index

242 { return m_zref[lev].get(); }

MOSTAverage::m_zref

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

Definition: ERF_MOSTAverage.H:348

Referenced by SurfaceLayer::get_zref().

Here is the caller graph for this function:

◆ make_MOSTAverage_at_level()

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

Make MOST-average data structures at one level.

Parameters

[in]	lev	level index
[in]	vars_old	old-time state and velocity fields
[in]	Theta_prim	primitive potential-temperature field
[in]	Qv_prim	primitive water-vapor field
[in]	Qr_prim	primitive rain-water field
[in]	z_phys_nd	nodal physical-height field

Make MOSTAverage data structures at one level.

Parameters

[in]	lev	Current level
[in]	vars_old	State data used by the average calculator
[in]	Theta_prim	Primitive theta component at this level
[in]	Qv_prim	Primitive water-vapor component at this level
[in]	Qr_prim	Primitive rain-water component at this level
[in]	z_phys_nd	Nodal physical height at this level

 {
     m_fields[lev].resize(m_nvar);
     m_rot_fields[lev].resize(m_nvar-1);
     m_averages[lev].resize(m_navg);
     m_z_phys_nd[lev] = z_phys_nd.get();
  
     bool use_terrain_fitted_coords = ( (m_terrain_type == TerrainType::StaticFittedMesh) ||
                                        (m_terrain_type == TerrainType::MovingFittedMesh) );
  
     bool use_eb = (m_terrain_type == TerrainType::EB);
  
     { // Nodal in x
         auto& mf = *vars_old[Vars::xvel];
         // Create a 2D ba, dm, & ghost cells
         const BoxArray& ba = mf.boxArray();
         BoxList bl2d = ba.boxList();
         for (auto& b : bl2d) { b.setRange(2,0); }
         BoxArray ba2d(std::move(bl2d));
         const DistributionMapping& dm = mf.DistributionMap();
         const int ncomp = 1;
         IntVect ng = mf.nGrowVect(); ng[2]=0;
  
         m_fields[lev][0] = vars_old[Vars::xvel];
         m_averages[lev][0] = std::make_unique<MultiFab>(ba2d,dm,ncomp,ng);
         m_averages[lev][0]->setVal(bogus_large_value);
         if (m_rotate) {
             m_rot_fields[lev][0] = std::make_unique<MultiFab>(ba,dm,ncomp,ng);
             MultiFab::Copy(*m_rot_fields[lev][0],mf,0,0,1,ng);
         } else {
             m_rot_fields[lev][0] = nullptr;
         }
     }
     { // Nodal in y
         auto& mf  = *vars_old[Vars::yvel];
         // Create a 2D ba, dm, & ghost cells
         const BoxArray& ba = mf.boxArray();
         BoxList bl2d = ba.boxList();
         for (auto& b : bl2d) { b.setRange(2,0); }
         BoxArray ba2d(std::move(bl2d));
         const DistributionMapping& dm = mf.DistributionMap();
         const int ncomp = 1;
         IntVect ng = mf.nGrowVect(); ng[2]=0;
  
         m_fields[lev][1] = vars_old[Vars::yvel];
         m_averages[lev][1] = std::make_unique<MultiFab>(ba2d,dm,ncomp,ng);
         m_averages[lev][1]->setVal(bogus_large_value);
         if (m_rotate) {
             m_rot_fields[lev][1] = std::make_unique<MultiFab>(ba,dm,ncomp,ng);
             MultiFab::Copy(*m_rot_fields[lev][1],mf,0,0,1,ng);
         } else {
             m_rot_fields[lev][1] = nullptr;
         }
     }
     { // CC vars
         auto& mf  = *Theta_prim;
         // Create a 2D ba, dm, & ghost cells
         const BoxArray& ba = mf.boxArray();
         BoxList bl2d = ba.boxList();
         for (auto& b : bl2d) { b.setRange(2,0); }
         BoxArray ba2d(std::move(bl2d));
         const DistributionMapping& dm = mf.DistributionMap();
         const int ncomp  = 1;
         const int incomp = 1;
         IntVect ng = mf.nGrowVect(); ng[2]=0;
  
         // Get field pointers
         m_fields[lev][2] = Theta_prim.get();
         m_fields[lev][3] = Qv_prim.get();
         m_fields[lev][4] = Qr_prim.get();
  
         // Initialize remaining multifabs
         for (int iavg(2); iavg < m_navg; ++iavg) {
             m_averages[lev][iavg] = std::make_unique<MultiFab>(ba2d,dm,ncomp,ng);
             m_averages[lev][iavg]->setVal(bogus_large_value);
         }
  
         // Default to dry
         m_averages[lev][3]->setVal(0.0);
  
         if (m_rotate) {
             m_rot_fields[lev][2] = std::make_unique<MultiFab>(ba,dm,ncomp,ng);
             m_rot_fields[lev][3] = std::make_unique<MultiFab>(ba,dm,ncomp,ng);
             MultiFab::Copy(*m_rot_fields[lev][2],*Theta_prim,0,0,1,ng);
             if (Qv_prim) MultiFab::Copy(*m_rot_fields[lev][3],*Qv_prim,0,0,1,ng);
         } else {
             m_rot_fields[lev][2] = nullptr;
             m_rot_fields[lev][3] = nullptr;
         }
  
         // Default zref to 10 and fill will true values later
         m_zref[lev] = std::make_unique<MultiFab>(ba2d,dm,1,ng);
         m_zref[lev]->setVal(zref_default);
  
         if (use_terrain_fitted_coords && m_norm_vec && m_interp) {
             m_x_pos[lev] = std::make_unique<MultiFab>(ba2d,dm,ncomp,ng);
             m_y_pos[lev] = std::make_unique<MultiFab>(ba2d,dm,ncomp,ng);
             m_z_pos[lev] = std::make_unique<MultiFab>(ba2d,dm,ncomp,ng);
         } else if (use_terrain_fitted_coords && m_interp) {
             m_x_pos[lev] = std::make_unique<MultiFab>(ba2d,dm,ncomp,ng);
             m_y_pos[lev] = std::make_unique<MultiFab>(ba2d,dm,ncomp,ng);
             m_z_pos[lev] = std::make_unique<MultiFab>(ba2d,dm,ncomp,ng);
         } else if (use_terrain_fitted_coords && m_norm_vec) {
             m_i_indx[lev] = std::make_unique<iMultiFab>(ba2d,dm,incomp,ng);
             m_j_indx[lev] = std::make_unique<iMultiFab>(ba2d,dm,incomp,ng);
             m_k_indx[lev] = std::make_unique<iMultiFab>(ba2d,dm,incomp,ng);
         } else {
             if (!use_eb) {
                 m_k_indx[lev] = std::make_unique<iMultiFab>(ba2d,dm,incomp,ng);
             }
         }
     }
     // Nodal in z (only used with terrain stress rotations)
     m_fields[lev][5] = vars_old[Vars::zvel];
  
     // Setup auxiliary data for spatial configuration & policy
     //--------------------------------------------------------
     if (use_terrain_fitted_coords && m_norm_vec && m_interp) { // Terrain w/ norm & w/ interpolation
         set_norm_positions_T(lev);
     } else if (use_terrain_fitted_coords && m_interp) {        // Terrain w/ interpolation
         set_z_positions_T(lev);
     } else if (use_terrain_fitted_coords && m_norm_vec) {      // Terrain w/ norm & w/o interpolation
         set_norm_indices_T(lev);
     } else if (use_terrain_fitted_coords) {                    // Terrain
         set_k_indices_T(lev);
     } else if (use_eb) {                                       // EB
         set_z_positions_EB(lev);
     } else {                                                   // No Terrain
         set_k_indices_N(lev);
     }
  
     // Setup normalization data for the chosen policy
     //--------------------------------------------------------
     switch(m_policy) {
     case 0: // Plane average
         set_plane_normalization(lev);
         break;
     case 1: // Local region/point
         set_region_normalization(lev);
         break;
     case 2: // EB
         set_eb_normalization(lev);
         break;
     default:
         AMREX_ALWAYS_ASSERT_WITH_MESSAGE(false, "Unknown policy for MOSTAverage!");
     }
  
     // Set up the exponential time filtering
     //--------------------------------------------------------
     if (m_t_avg) {
         // Exponential filter function
         m_fact_old = static_cast<amrex::Real>(std::exp(-1.0 / m_time_window));
  
         // Enforce discrete normalization: (mfn*val_new + mfo*val_old)
         m_fact_new = one - m_fact_old;
  
         // None of the averages are initialized
         m_t_init.resize(m_maxlev,0);
     }
  
     // Correction to the mean surface velocity at this level
     m_Vsg[lev] = zero;
     if (include_subgrid_vel) {
         Print() << "Subgrid velocity scale correction at level : " << lev << ' ';
         const auto dxArr = m_geom[lev].CellSizeArray();
         Real dx = std::sqrt(dxArr[0]*dxArr[1]);
         if (dx > Real(5000.)) {
             m_Vsg[lev] = Real(0.32) * std::pow(dx/Real(5000.)-1, Real(0.33));
         }
         Print() << m_Vsg[lev] << std::endl;
     }
 }

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] other source object

◆ operator=() [2/2]

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

deletenoexcept

Deleted move-assignment operator.

Parameters

[in] other source object

◆ set_eb_normalization()

void MOSTAverage::set_eb_normalization ( const int & lev )

Compute total embedded-boundary surface area.

Parameters

[in] lev level index

Function to compute normalization for average over EB.

Parameters

[in] lev Current level

 {
     AMREX_ALWAYS_ASSERT(m_eb_vec[lev] != nullptr);
  
     // Get EB data - need both cell-centered and face-centered flags
     const auto& cc_flags = m_eb_vec[lev]->get_const_factory()->getMultiEBCellFlagFab();
  
     // Get area fractions for different centerings
     auto cc_afrac = m_eb_vec[lev]->get_const_factory()->getAreaFrac();  // Cell-centered area fractions (Array of 3 MultiCutFab*)
  
     // Initialize storage
     m_total_bndry_area.resize(m_maxlev);
     m_total_bndry_area[lev].resize(m_navg, zero);
     m_plane_average.resize(m_maxlev);
     m_plane_average[lev].resize(m_navg, zero);
  
     // Compute total area for each field type based on its centering
     // iavg: 0=U(xface), 1=V(yface), 2=T(cc), 3=Qv(cc), 4=Tv(cc), 5=Umag(cc)
  
     // Get geometry for cell sizes
     auto const& dx_arr = m_geom[lev].CellSizeArray();
     Real dx = dx_arr[0];
     Real dy = dx_arr[1];
     Real dz = dx_arr[2];
  
     // GPU array to accumulate areas
     Gpu::DeviceVector<Real> area_vec(m_navg, zero);
     Real* area_device = area_vec.data();
  
     // All fields are now averaged on cell-centered grid
     // Compute total area once on cell-centered grid and use for all iavg
     Real total_area = zero;
 #ifdef _OPENMP
 #pragma omp parallel if (Gpu::notInLaunchRegion())
 #endif
     for (MFIter mfi(cc_flags, TileNoZ()); mfi.isValid(); ++mfi) {
         const auto& flag = cc_flags[mfi];
  
         // Skip boxes that are not singlevalued (MultiCutFab only has data for singlevalued boxes)
         if (flag.getType() != FabType::singlevalued) continue;
  
         Box bx = mfi.tilebox();
         auto const flag_arr = flag.const_array();
         auto const afrac_x = cc_afrac[0]->const_array(mfi);
         auto const afrac_y = cc_afrac[1]->const_array(mfi);
         auto const afrac_z = cc_afrac[2]->const_array(mfi);
  
         // Sum area only for cut cells using atomic reduction
         ParallelFor(Gpu::KernelInfo().setReduction(true), bx, [=]
         AMREX_GPU_DEVICE(int i, int j, int k, Gpu::Handler const& handler) noexcept
         {
             if (flag_arr(i,j,k).isSingleValued()) {
                 // Compute area from face-centered area fractions
                 Real axm = afrac_x(i  ,j  ,k  );
                 Real axp = afrac_x(i+1,j  ,k  );
                 Real aym = afrac_y(i  ,j  ,k  );
                 Real ayp = afrac_y(i  ,j+1,k  );
                 Real azm = afrac_z(i  ,j  ,k  );
                 Real azp = afrac_z(i  ,j  ,k+1);
  
                 Real adx = (axm - axp) * dy * dz;
                 Real ady = (aym - ayp) * dx * dz;
                 Real adz = (azm - azp) * dx * dy;
  
                 Real area = std::sqrt(adx*adx + ady*ady + adz*adz);
  
                 Gpu::deviceReduceSum(&area_device[0], area, handler);
             }
         });
     }
  
     // Copy to host and sum across MPI ranks
     Gpu::copy(Gpu::deviceToHost, area_vec.begin(), area_vec.begin() + 1, &total_area);
     ParallelDescriptor::ReduceRealSum(&total_area, 1);
  
     // Set the same total area for all iavg
     for (int iavg = 0; iavg < m_navg; ++iavg) {
         m_total_bndry_area[lev][iavg] = total_area;
     }
  
     Print() << "EB surface area on cell-centerd grid at level " << lev << ": " << total_area << std::endl;
 }

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] lev level index

Function to set K indices without terrain.

Parameters

[in] lev Current level

 {
     ParmParse pp(m_pp_prefix);
     Real zref_tmp = zref_default;
     auto read_z = pp.query("most.zref",zref_tmp);
     auto read_k = pp.queryarr("most.k_arr_in",m_k_in);
  
     // Default behavior is to use the first cell center
     if (!read_z && !read_k) {
         Real m_zlo = m_geom[0].ProbLo(2);
         Real m_dz  = m_geom[0].CellSize(2);
         zref_tmp = m_zlo + myhalf * m_dz;
         m_zref[lev]->setVal( zref_tmp );
         Print() << "Reference height for MOST set to " << zref_tmp << std::endl;
         read_z = true;
     }
  
     // Specify z_ref & compute k_indx (z_ref takes precedence)
     if (read_z) {
         Real m_zlo = m_geom[lev].ProbLo(2);
         Real m_zhi = m_geom[lev].ProbHi(2);
         Real m_dz  = m_geom[lev].CellSize(2);
  
         amrex::ignore_unused(m_zhi);
  
         AMREX_ALWAYS_ASSERT_WITH_MESSAGE(zref_tmp >= m_zlo + myhalf * m_dz,
                                          "Query point must be past first z-cell!");
  
         AMREX_ALWAYS_ASSERT_WITH_MESSAGE(zref_tmp <= m_zhi - myhalf * m_dz,
                                          "Query point must be below the last z-cell!");
  
         int lk = static_cast<int>(floor((zref_tmp - m_zlo) / m_dz - myhalf));
  
         m_zref[lev]->setVal( (lk + myhalf) * m_dz + m_zlo );
  
         AMREX_ALWAYS_ASSERT(lk >= m_radius);
  
         m_k_indx[lev]->setVal(lk);
     // Specified k_indx & compute z_ref
     } else if (read_k) {
         AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_k_in[lev] >= m_radius,
                                   "K index must be larger than averaging radius!");
         m_k_indx[lev]->setVal(m_k_in[lev]);
  
         // TODO: check that z_ref is constant across levels
         Real m_zlo = m_geom[0].ProbLo(2);
         Real m_dz  = m_geom[0].CellSize(2);
         m_zref[lev]->setVal( ((Real)m_k_in[0] + myhalf) * m_dz + m_zlo );
     }
 }

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 )

Populate the 2D k-index iMultiFab with terrain.

Parameters

[in] lev level index

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

Parameters

[in] lev Current level

 {
     // Peel back the level
     auto& fields = m_fields[lev];
  
     // MFIter over CC data
     int imf_cc = 2;
  
     ParmParse pp(m_pp_prefix);
     Real zref_tmp = zref_default;
     auto read_z = pp.query("most.zref",zref_tmp);
     auto read_k = pp.queryarr("most.k_arr_in",m_k_in);
     int klo     = m_geom[lev].Domain().smallEnd(2);
  
     // Allow default zref
     if (!read_z) {
         Print() << "most.zref not specified, query distance default is " << zref_tmp << std::endl;
         read_z = true;
     }
  
     // No default behavior with terrain (we can't tell the difference between
     // vertical grid stretching and true terrain)
     AMREX_ALWAYS_ASSERT_WITH_MESSAGE(read_z != read_k,
                                      "Need to specify zref or k_arr_in for MOST");
  
     // Capture for device
     Real d_zref   = zref_tmp;
     Real d_radius = static_cast<Real>(m_radius);
     amrex::ignore_unused(d_radius);
  
     // Specify z_ref & compute k_indx (z_ref takes precedence)
     if (read_z) {
         int kmax = m_geom[lev].Domain().bigEnd(2);
         for (MFIter mfi(*fields[imf_cc], TileNoZ()); mfi.isValid(); ++mfi) {
             Box npbx = mfi.tilebox(IntVect(1,1,0),IntVect(1,1,0));
  
             if (npbx.smallEnd(2) != klo) { continue; }
  
             npbx.makeSlab(2,klo);
  
             const auto z_phys_arr = m_z_phys_nd[lev]->const_array(mfi);
             auto k_arr = m_k_indx[lev]->array(mfi);
             auto zref_arr = m_zref[lev]->array(mfi);
             ParallelFor(npbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
             {
                 k_arr(i,j,k) = klo;
                 bool found = false;
                 Real z_bot_face  = fourth * ( z_phys_arr(i  ,j  ,k) + z_phys_arr(i+1,j  ,k)
                                           + z_phys_arr(i  ,j+1,k) + z_phys_arr(i+1,j+1,k) );
                 Real z_target    = z_bot_face + d_zref;
                 for (int lk(klo); lk<=kmax; ++lk) {
                     Real z_lo = fourth * ( z_phys_arr(i,j  ,lk  ) + z_phys_arr(i+1,j  ,lk  )
                                        + z_phys_arr(i,j+1,lk  ) + z_phys_arr(i+1,j+1,lk  ) );
                     Real z_hi = fourth * ( z_phys_arr(i,j  ,lk+1) + z_phys_arr(i+1,j  ,lk+1)
                                        + z_phys_arr(i,j+1,lk+1) + z_phys_arr(i+1,j+1,lk+1) );
                     if (z_target > z_lo && z_target < z_hi){
                         AMREX_ALWAYS_ASSERT_WITH_MESSAGE(lk >= d_radius,
                                                          "K index must be larger than averaging radius!");
                         k_arr(i,j,0) = lk;
                         zref_arr(i,j,0) = myhalf * (z_hi + z_lo) - z_bot_face;
                         found = true;
                         break;
                     }
                 }
                 AMREX_ALWAYS_ASSERT_WITH_MESSAGE(found,
                                                  "zref not found with terrain!");
             });
         }
     // Specified k_indx & compute z_ref
     } else if (read_k) {
         AMREX_ALWAYS_ASSERT_WITH_MESSAGE(false, "Specified k-indx with terrain not implemented!");
     }
 }

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] lev level index

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

Parameters

[in] lev Current level

 {
     // Peel back the level
     auto& fields = m_fields[lev];
  
     // MFIter over CC data
     int imf_cc = 2;
  
     ParmParse pp(m_pp_prefix);
     Real zref_tmp = zref_default;
     auto read_zref = pp.query("most.zref",zref_tmp);
     if (!read_zref) {
         Print() << "most.zref not specified, query distance default is " << zref_tmp << std::endl;
     }
     int klo = m_geom[lev].Domain().smallEnd(2);
  
     // Capture for device
     Real d_zref   = zref_tmp;
     Real d_radius = static_cast<Real>(m_radius);
  
     const auto dxInv  = m_geom[lev].InvCellSizeArray();
     IntVect ng = m_k_indx[lev]->nGrowVect(); ng[2]=0;
     for (MFIter mfi(*fields[imf_cc], TileNoZ()); mfi.isValid(); ++mfi) {
         Box npbx  = mfi.tilebox(IntVect(1,1,0),IntVect(1,1,0));
  
         if (npbx.smallEnd(2) != klo) { continue; }
  
         int kmax = npbx.bigEnd(2);
  
         npbx.makeSlab(2,klo);
  
         const auto z_phys_arr = m_z_phys_nd[lev]->const_array(mfi);
         auto i_arr = m_i_indx[lev]->array(mfi);
         auto j_arr = m_j_indx[lev]->array(mfi);
         auto k_arr = m_k_indx[lev]->array(mfi);
         auto zref_arr = m_zref[lev]->array(mfi);
         ParallelFor(npbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
         {
             // Elements of normal vector
             Real met_h_xi  = Compute_h_xi_AtCellCenter (i,j,k,dxInv,z_phys_arr);
             Real met_h_eta = Compute_h_eta_AtCellCenter(i,j,k,dxInv,z_phys_arr);
             Real mag = std::sqrt(met_h_xi*met_h_xi + met_h_eta*met_h_eta + one);
  
             // Unit-normal vector scaled by z_ref
             Real delta_x = -met_h_xi/mag  * d_zref;
             Real delta_y = -met_h_eta/mag * d_zref;
             Real delta_z = one/mag * d_zref;
  
             // Compute i & j as displacements (no grid stretching)
             int delta_i  = static_cast<int>(std::round(delta_x*dxInv[0]));
             int delta_j  = static_cast<int>(std::round(delta_y*dxInv[1]));
             int i_new    = i + delta_i;
             int j_new    = j + delta_j;
             i_arr(i,j,0) = i_new;
             j_arr(i,j,0) = j_new;
  
             // Search for k (grid is stretched in z)
             Real z_bot_face  = fourth * ( z_phys_arr(i  ,j  ,k) + z_phys_arr(i+1,j  ,k)
                                           + z_phys_arr(i  ,j+1,k) + z_phys_arr(i+1,j+1,k) );
             Real z_target    = z_bot_face + delta_z;
             k_arr(i,j,0)     = klo;
             zref_arr(i,j,0)  = myhalf * z_bot_face +
                                Real(0.125) * ( z_phys_arr(i  ,j  ,k+1) + z_phys_arr(i+1,j  ,k+1)
                                              + z_phys_arr(i  ,j+1,k+1) + z_phys_arr(i+1,j+1,k+1) );
             for (int lk(klo); lk<=kmax; ++lk) {
                 Real z_lo = fourth * ( z_phys_arr(i_new,j_new  ,lk  ) + z_phys_arr(i_new+1,j_new  ,lk  )
                                        + z_phys_arr(i_new,j_new+1,lk  ) + z_phys_arr(i_new+1,j_new+1,lk  ) );
                 Real z_hi = fourth * ( z_phys_arr(i_new,j_new  ,lk+1) + z_phys_arr(i_new+1,j_new  ,lk+1)
                                        + z_phys_arr(i_new,j_new+1,lk+1) + z_phys_arr(i_new+1,j_new+1,lk+1) );
                 if (z_target > z_lo && z_target < z_hi){
                     AMREX_ALWAYS_ASSERT_WITH_MESSAGE(lk >= d_radius,
                                                      "K index must be larger than averaging radius!");
                     amrex::ignore_unused(d_radius);
                     k_arr(i,j,0) = lk;
                     zref_arr(i,j,0) = myhalf * (z_hi + z_lo) - z_bot_face;
                     break;
                 }
             }
         });
     }
 }

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] lev level index

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

Parameters

[in] lev Current level

 {
     // Peel back the level
     auto& fields = m_fields[lev];
  
     // MFIter over CC data
     int imf_cc = 2;
  
     ParmParse pp(m_pp_prefix);
     Real zref_tmp = zref_default;
     auto read_zref = pp.query("most.zref",zref_tmp);
     if (!read_zref) {
         Print() << "most.zref not specified, query distance default is " << zref_tmp << std::endl;
     }
     int klo = m_geom[lev].Domain().smallEnd(2);
  
     // Capture for device
     Real d_zref = zref_tmp;
     const auto plo = m_geom[lev].ProbLoArray();
  
     RealVect base;
     const auto dx = m_geom[lev].CellSizeArray();
     const auto dxInv  = m_geom[lev].InvCellSizeArray();
     IntVect ng = m_x_pos[lev]->nGrowVect(); ng[2]=0;
     const int position_ng = (m_radius > 1) ? m_radius : 1;
     for (MFIter mfi(*fields[imf_cc], TileNoZ()); mfi.isValid(); ++mfi) {
         Box npbx  = mfi.tilebox(IntVect(1,1,0),IntVect(position_ng,position_ng,0));
         Box gtbx  = mfi.growntilebox(ng);
         RealBox grb{gtbx,dx.data(),base.dataPtr()};
  
         if (npbx.smallEnd(2) != klo) { continue; }
  
         npbx.makeSlab(2,klo);
  
         const auto z_phys_arr = m_z_phys_nd[lev]->const_array(mfi);
         auto x_pos_arr   = m_x_pos[lev]->array(mfi);
         auto y_pos_arr   = m_y_pos[lev]->array(mfi);
         auto z_pos_arr   = m_z_pos[lev]->array(mfi);
         auto zref_arr    = m_zref[lev]->array(mfi);
         ParallelFor(npbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
         {
             // Elements of normal vector
             Real met_h_xi  = Compute_h_xi_AtCellCenter (i,j,k,dxInv,z_phys_arr);
             Real met_h_eta = Compute_h_eta_AtCellCenter(i,j,k,dxInv,z_phys_arr);
             Real imag = one / std::sqrt(met_h_xi*met_h_xi + met_h_eta*met_h_eta + one);
  
             // Unit-normal vector scaled by z_ref
             Real delta_x = -met_h_xi  * imag * d_zref;
             Real delta_y = -met_h_eta * imag * d_zref;
             Real delta_z =              imag * d_zref;
  
             // Position of the current node (indx:0,0,1)
             Real x0 = plo[0] + ((Real) i + myhalf) * dx[0];
             Real y0 = plo[1] + ((Real) j + myhalf) * dx[1];
  
             // Final position at end of vector
             x_pos_arr(i,j,0) = x0 + delta_x;
             y_pos_arr(i,j,0) = y0 + delta_y;
             Real z_bot_face  = fourth * ( z_phys_arr(i  ,j  ,k) + z_phys_arr(i+1,j  ,k)
                                           + z_phys_arr(i  ,j+1,k) + z_phys_arr(i+1,j+1,k) );
             z_pos_arr(i,j,0) = z_bot_face + delta_z;
  
             // NOTE: Normal vector end point can be below the surface for concave regions.
             //       Here we protect against that by augmenting the normal if needed.
             int i_new = (int) ((x_pos_arr(i,j,0) - plo[0]) / dx[0] - myhalf);
             int j_new = (int) ((y_pos_arr(i,j,0) - plo[1]) / dx[1] - myhalf);
             Real z_new_bot_face = fourth * ( z_phys_arr(i_new,j_new  ,k) + z_phys_arr(i_new+1,j_new  ,k)
                                              + z_phys_arr(i_new,j_new+1,k) + z_phys_arr(i_new+1,j_new+1,k) );
             if (z_pos_arr(i,j,0) < z_new_bot_face) {
                 z_pos_arr(i,j,0) = z_new_bot_face + delta_z;
             }
  
             zref_arr(i,j,0) = delta_z;
  
             // Destination position must be contained on the current process!
             Real pos[] = {x_pos_arr(i,j,0)-plo[0],y_pos_arr(i,j,0)-plo[1],myhalf*dx[2]};
             amrex::ignore_unused(pos);
             AMREX_ALWAYS_ASSERT_WITH_MESSAGE(grb.contains(&pos[0]),
                                               "Query point outside of proc domain!");
         });
     }
 }

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 )

Compute number of cells per averaging plane.

Parameters

[in] lev level index

Function to compute normalization for plane average.

Parameters

[in] lev Current level

 {
     // Cells per plane and temp avg storage
     m_ncell_plane.resize(m_maxlev);
     m_plane_average.resize(m_maxlev);
  
     // True domain not used for normalization
     Box domain = m_geom[lev].Domain();
  
     // NOTE: Level 0 spans the whole domain, but finer
     //       levels do not have such a restriction.
     //       For now, use the bounding box of the boxArray
     //       for normalization, consistent with avg routine.
  
     // Bounded box of CC data used for normalization
     Box bnd_bx = (m_fields[lev][2]->boxArray()).minimalBox();
  
     // NOTE: Bounding box must lie on the periodic boundaries
     //       in order to trip the is_per flag
  
     // Num components, plane avg, cells per plane
     Array<int,AMREX_SPACEDIM> is_per = {0,0,0};
     for (int idim(0); idim < AMREX_SPACEDIM-1; ++idim) {
         if ( m_geom[lev].isPeriodic(idim) &&
              bnd_bx.bigEnd(idim)==domain.bigEnd(idim) &&
              bnd_bx.smallEnd(idim)==domain.smallEnd(idim) ) { is_per[idim] = 1; }
     }
  
     m_ncell_plane[lev].resize(m_navg);
     m_plane_average[lev].resize(m_navg);
     for (int iavg(0); iavg < m_navg; ++iavg) {
         // Convert bnd_bx to current index type
         IndexType ixt = m_averages[lev][iavg]->boxArray().ixType();
         bnd_bx.convert(ixt);
         IntVect bnd_bx_lo(bnd_bx.loVect());
         IntVect bnd_bx_hi(bnd_bx.hiVect());
  
         m_plane_average[lev][iavg] = zero;
  
         m_ncell_plane[lev][iavg] = 1;
         for (int idim(0); idim < AMREX_SPACEDIM; ++idim) {
             if (idim != 2) {
                 if (ixt.nodeCentered(idim) && is_per[idim]) {
                     m_ncell_plane[lev][iavg] *= (bnd_bx_hi[idim] - bnd_bx_lo[idim]);
                 } else {
                     m_ncell_plane[lev][iavg] *= (bnd_bx_hi[idim] - bnd_bx_lo[idim] + 1);
                 }
             }
         } // idim
     } // iavg
 }

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] lev level index

127 {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] lev level index

Function to set the rotated velocities.

Parameters

[in] lev Current level

 {
     // Peel back the level
     auto& fields     = m_fields[lev];
     auto& rot_fields = m_rot_fields[lev];
     auto  z_phys_nd  = m_z_phys_nd[lev];
  
     // Inverse grid size
     const auto dxInv = m_geom[lev].InvCellSizeArray();
  
     // Single MFIter over CC data
     int imf_cc = 2;
  
     // Populate rotated U & V for terrain
 #ifdef _OPENMP
 #pragma omp parallel if (Gpu::notInLaunchRegion())
 #endif
     for (MFIter mfi(*fields[imf_cc], TileNoZ()); mfi.isValid(); ++mfi) {
         Box ubx = mfi.tilebox(IntVect(1,0,0));
         Box vbx = mfi.tilebox(IntVect(0,1,0));
  
         const Array4<const Real>& z_phys_arr = z_phys_nd->const_array(mfi);
  
         const Array4<const Real>& u_arr = fields[0]->const_array(mfi);
         const Array4<const Real>& v_arr = fields[1]->const_array(mfi);
         const Array4<const Real>& w_arr = fields[5]->const_array(mfi);
  
         const Array4<Real>& u_rot_arr = rot_fields[0]->array(mfi);
         const Array4<Real>& v_rot_arr = rot_fields[1]->array(mfi);
  
         // U rotated magnitude
         ParallelFor(ubx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
         {
             // Elements of first tangent vector
             Real met_h_xi    = Compute_h_xi_AtIface(i,j,k,dxInv,z_phys_arr);
             u_rot_arr(i,j,k) = (u_arr(i,j,k) + met_h_xi*w_arr(i,j,k))
                              / std::sqrt(met_h_xi*met_h_xi + one);
         });
  
         // V rotated magnitude
         ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
         {
             // Elements of second tangent vector
             Real met_h_eta   = Compute_h_eta_AtJface(i,j,k,dxInv,z_phys_arr);
             v_rot_arr(i,j,k) = (v_arr(i,j,k) + met_h_eta*w_arr(i,j,k))
                              / std::sqrt(met_h_eta*met_h_eta + one);
         });
     }
  
     // Direct copy of other scalar variables
     MultiFab::Copy(*rot_fields[2],*fields[2],0,0,1,rot_fields[2]->nGrowVect());
     if (fields[3]) MultiFab::Copy(*rot_fields[3],*fields[3],0,0,1,rot_fields[3]->nGrowVect());
 }

Referenced by compute_averages().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ set_z_positions_EB()

void MOSTAverage::set_z_positions_EB ( const int & lev )

Populate embedded-boundary z positions.

Parameters

[in] lev level index

Function to set K indices for EB.

Parameters

[in] lev Current level

 {
     Real zref_tmp = zref_default;
     ParmParse pp(m_pp_prefix);
     auto read_z = pp.query("most.zref",zref_tmp);
  
     if (read_z) {
         m_zref[lev]->setVal( zref_tmp );
     // Default behavior is to use the first cell center
     } else {
         Real m_dz  = m_geom[0].CellSize(2);
         zref_tmp = myhalf * m_dz;
         m_zref[lev]->setVal( zref_tmp );
         Print() << "Reference height for MOST set to " << zref_tmp << std::endl;
     }
 }

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] lev level index

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

Parameters

[in] lev Current level

 {
     // Peel back the level
     auto& fields = m_fields[lev];
  
     // MFIter over CC data
     int imf_cc = 2;
  
     ParmParse pp(m_pp_prefix);
     Real zref_tmp = zref_default;
     auto read_zref = pp.query("most.zref",zref_tmp);
     if (!read_zref) {
         Print() << "most.zref not specified, query distance default is " << zref_tmp << std::endl;
     } else {
         m_zref[lev]->setVal(zref_tmp);
     }
     int klo = m_geom[lev].Domain().smallEnd(2);
  
     // Capture for device
     Real d_zref = zref_tmp;
     const auto plo = m_geom[lev].ProbLoArray();
  
     RealVect base;
     const auto dx = m_geom[lev].CellSizeArray();
     IntVect ng = m_x_pos[lev]->nGrowVect(); ng[2]=0;
     const int position_ng = (m_radius > 1) ? m_radius : 1;
     for (MFIter mfi(*fields[imf_cc], TileNoZ()); mfi.isValid(); ++mfi) {
         Box npbx  = mfi.tilebox(IntVect(1,1,0),IntVect(position_ng,position_ng,0));
         Box gtbx  = mfi.growntilebox(ng);
  
         if (npbx.smallEnd(2) != klo) { continue; }
  
         npbx.makeSlab(2,klo);
  
         RealBox grb{gtbx,dx.data(),base.dataPtr()};
  
         const auto z_phys_arr = m_z_phys_nd[lev]->const_array(mfi);
         auto x_pos_arr = m_x_pos[lev]->array(mfi);
         auto y_pos_arr = m_y_pos[lev]->array(mfi);
         auto z_pos_arr = m_z_pos[lev]->array(mfi);
         ParallelFor(npbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
         {
             // Final position at end of vector
             x_pos_arr(i,j,0) = plo[0] + ((Real) i + myhalf) * dx[0];
             y_pos_arr(i,j,0) = plo[1] + ((Real) j + myhalf) * dx[1];
             Real z_bot_face  = fourth * ( z_phys_arr(i  ,j  ,k) + z_phys_arr(i+1,j  ,k)
                                           + z_phys_arr(i  ,j+1,k) + z_phys_arr(i+1,j+1,k) );
             z_pos_arr(i,j,0) = z_bot_face + d_zref;
  
             // Destination position must be contained on the current process!
             Real pos[] = {x_pos_arr(i,j,0)-plo[0],y_pos_arr(i,j,0)-plo[1],myhalf*dx[2]};
             amrex::ignore_unused(pos);
             AMREX_ALWAYS_ASSERT_WITH_MESSAGE(grb.contains(&pos[0]),
                                              "Query point outside of proc domain!");
         });
     }
 }

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]	xp	X-position
[in]	yp	Y-position
[in]	zp	Z-position
[out]	interp_vals	Values interpolated
[in]	interp_array	Array to interpolate on
[in]	z_arr	Physical heights
[in]	plo	Problem lower bounds
[in]	dxi	Inverse cell size array
[in]	interp_comp	Number of components to interpolate

     {
         // Search to get z/k
         bool found = false;
         int kmax   = ubound(z_arr).z;
         amrex::Real zval     = zero;
         amrex::Real z_target = zp;
  
         // Map position to i,j (must be same mapping in cpp file)
         amrex::Real ireal = (xp - plo[0]) * dxi[0];
         amrex::Real jreal = (yp - plo[1]) * dxi[1];
         int i_new = (int) (ireal - myhalf);
         int j_new = (int) (jreal - myhalf);
  
         for (int lk(0); lk<kmax; ++lk) {
             amrex::Real z_lo = fourth * ( z_arr(i_new,j_new  ,lk  ) + z_arr(i_new+1,j_new  ,lk  )
                                       + z_arr(i_new,j_new+1,lk  ) + z_arr(i_new+1,j_new+1,lk  ) );
             amrex::Real z_hi = fourth * ( z_arr(i_new,j_new  ,lk+1) + z_arr(i_new+1,j_new  ,lk+1)
                                       + z_arr(i_new,j_new+1,lk+1) + z_arr(i_new+1,j_new+1,lk+1) );
             if (z_target > z_lo && z_target < z_hi){
                 found = true;
                 zval  = (amrex::Real) lk + ((z_target - z_lo) / (z_hi - z_lo)) + myhalf;
                 break;
             }
         }
  
         amrex::ignore_unused(found);
         AMREX_ALWAYS_ASSERT_WITH_MESSAGE(found, "MOSTAverage: Height above terrain not found, try increasing z_ref!");
  
         // NOTE: This is the point ahead of the current i,j (e.g. i/j_new + 1)
         const amrex::RealVect lx(ireal + myhalf, jreal + myhalf, zval);
  
         const amrex::IntVect ijk = lx.floor();
  
         int i = ijk[0]; int j = ijk[1]; int k = ijk[2];
  
         // Convert ijk (IntVect) to a RealVect explicitly
         amrex::RealVect ijk_r(static_cast<amrex::Real>(ijk[0]),
                               static_cast<amrex::Real>(ijk[1]),
                               static_cast<amrex::Real>(ijk[2]));
  
         // Weights
         const amrex::RealVect sx_hi = lx  - ijk_r;
         const amrex::RealVect sx_lo = one - sx_hi;
  
         for (int n = 0; n < interp_comp; n++) {
             interp_vals[n] = sx_lo[0]*sx_lo[1]*sx_lo[2]*interp_array(i-1, j-1, k-1,n) +
                              sx_lo[0]*sx_lo[1]*sx_hi[2]*interp_array(i-1, j-1, k  ,n) +
                              sx_lo[0]*sx_hi[1]*sx_lo[2]*interp_array(i-1, j  , k-1,n) +
                              sx_lo[0]*sx_hi[1]*sx_hi[2]*interp_array(i-1, j  , k  ,n) +
                              sx_hi[0]*sx_lo[1]*sx_lo[2]*interp_array(i  , j-1, k-1,n) +
                              sx_hi[0]*sx_lo[1]*sx_hi[2]*interp_array(i  , j-1, k  ,n) +
                              sx_hi[0]*sx_hi[1]*sx_lo[2]*interp_array(i  , j  , k-1,n) +
                              sx_hi[0]*sx_hi[1]*sx_hi[2]*interp_array(i  , j  , k  ,n);
         }
     }

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]	lev	level index
[in]	vars_old	old-time state and velocity fields
[in]	Theta_prim	primitive potential-temperature fields by level
[in]	Qv_prim	primitive water-vapor fields by level
[in]	Qr_prim	primitive rain-water fields by level

Function to reset the pointers to field variables.

Parameters

[in]	lev	Current level
[in]	vars_old	Conserved variables at each level
[in]	Theta_prim	Primitive theta component at each level
[in]	Qv_prim	Primitive water-vapor component at each level
[in]	Qr_prim	Primitive rain-water component at each level

 {
     m_fields[lev][0] = &vars_old[lev][Vars::xvel];
     m_fields[lev][1] = &vars_old[lev][Vars::yvel];
     m_fields[lev][2] = Theta_prim[lev].get();
     m_fields[lev][3] = Qv_prim[lev].get();
     m_fields[lev][4] = Qr_prim[lev].get();
     m_fields[lev][5] = &vars_old[lev][Vars::zvel];
 }

Referenced by SurfaceLayer::update_mac_ptrs().

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] lev level index

Function to write averages to text file.

Parameters

[in] lev Current level

 {
     // Peel back the level
     auto& fields   = m_fields[lev];
     auto& averages = m_averages[lev];
  
     // MFIter on CC
     int imf_cc = 2;
  
     int klo = m_geom[lev].Domain().smallEnd(2);
  
     int navg = m_navg - 1;
  
     std::ofstream ofile;
     ofile.open ("MOST_averages.txt");
     ofile << "Averages computed via MOSTAverages class:\n";
  
     for (MFIter mfi(*fields[imf_cc], TileNoZ()); mfi.isValid(); ++mfi) {
         Box pbx = mfi.tilebox();
  
         if(pbx.smallEnd(2) != klo) { continue; }
  
         pbx.makeSlab(2,klo);
  
         int il = pbx.smallEnd(0); int iu = pbx.bigEnd(0);
         int jl = pbx.smallEnd(1); int ju = pbx.bigEnd(1);
  
         for (int j(jl); j <= ju; ++j) {
             for (int i(il); i <= iu; ++i) {
                 ofile << "(I,J): " << "(" << i << "," << j << ")" << "\n";
                 int k = 0;
                 for (int iavg(0); iavg <= navg; ++iavg) {
                     auto mf_arr = averages[iavg]->array(mfi);
                     ofile << "iavg val: "
                           << iavg << ' '
                           << mf_arr(i,j,k) << "\n";
                 }
                 ofile << "\n";
             }
         }
     }
     ofile.close();
 }

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] lev level index

Function to write the K indices to text file.

Parameters

[in] lev Current level

 {
     // Peel back the level
     auto& fields   = m_fields[lev];
     auto& k_indx   = m_k_indx[lev];
  
     // MFIter on CC
     int imf_cc = 2;
  
     int klo = m_geom[lev].Domain().smallEnd(2);
  
     std::ofstream ofile;
     ofile.open ("MOST_k_indices.txt");
     ofile << "K indices used to compute averages via MOSTAverages class:\n";
  
     for (MFIter mfi(*fields[imf_cc], TileNoZ()); mfi.isValid(); ++mfi) {
         Box pbx = mfi.tilebox();
  
         if(pbx.smallEnd(2) != klo) { continue; }
  
         pbx.makeSlab(2,klo);
  
         int il = pbx.smallEnd(0); int iu = pbx.bigEnd(0);
         int jl = pbx.smallEnd(1); int ju = pbx.bigEnd(1);
  
         auto k_arr = k_indx->array(mfi);
  
         for (int j(jl); j <= ju; ++j) {
             for (int i(il); i <= iu; ++i) {
                 ofile << "(I,J): " << "(" << i << "," << j << ")" << "\n";
                 int k = 0;
                 ofile << "K_ind: "
                       << k_arr(i,j,k) << "\n";
                 ofile << "\n";
             }
         }
     }
     ofile.close();
 }

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] lev level index

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

Parameters

[in] lev Current level

 {
     // Peel back the level
     auto& fields   = m_fields[lev];
     auto& k_indx   = m_k_indx[lev];
     auto& j_indx   = m_j_indx[lev];
     auto& i_indx   = m_i_indx[lev];
  
     // MFIter on CC
     int imf_cc = 2;
  
     int klo = m_geom[lev].Domain().smallEnd(2);
  
     std::ofstream ofile;
     ofile.open ("MOST_ijk_indices.txt");
     ofile << "IJK indices used to compute averages via MOSTAverages class:\n";
  
     for (MFIter mfi(*fields[imf_cc], TileNoZ()); mfi.isValid(); ++mfi) {
         Box pbx = mfi.tilebox();
  
         if(pbx.smallEnd(2) != klo) { continue; }
  
         pbx.makeSlab(2,klo);
  
         int il = pbx.smallEnd(0); int iu = pbx.bigEnd(0);
         int jl = pbx.smallEnd(1); int ju = pbx.bigEnd(1);
  
         auto k_arr = k_indx->array(mfi);
         auto j_arr = j_indx ? j_indx->array(mfi) : Array4<int> {};
         auto i_arr = i_indx ? i_indx->array(mfi) : Array4<int> {};
  
         for (int j(jl); j <= ju; ++j) {
             for (int i(il); i <= iu; ++i) {
                 ofile << "(I1,J1,K1): " << "(" << i << "," << j << "," << 0 << ")" << "\n";
  
                 int k  = 0;
                 int km = k_arr(i,j,k);
                 int jm = j_arr ? j_arr(i,j,k) : j;
                 int im = i_arr ? i_arr(i,j,k) : i;
  
                 ofile << "(I2,J2,K2): "
                       << "(" << im << "," << jm << "," << km << ")" << "\n";
                 ofile << "\n";
             }
         }
     }
     ofile.close();
 }

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]	lev	level index
[in]	j	y-index of the plane to write

Function to write X & Z positions to text file.

Parameters

[in]	lev	Current level
[in]	j	Index in y-dir

 {
     // Peel back the level
     auto& fields   = m_fields[lev];
     auto& x_pos_mf = m_x_pos[lev];
     auto& z_pos_mf = m_z_pos[lev];
  
     // MFIter on CC
     int imf_cc = 2;
  
     int klo = m_geom[lev].Domain().smallEnd(2);
  
     std::ofstream ofile;
     ofile.open ("MOST_xz_positions.txt");
  
     for (MFIter mfi(*fields[imf_cc], TileNoZ()); mfi.isValid(); ++mfi) {
         Box pbx = mfi.tilebox();
  
         if(pbx.smallEnd(2) != klo) { continue; }
  
         pbx.makeSlab(2,klo);
  
         int il = pbx.smallEnd(0); int iu = pbx.bigEnd(0);
  
         auto x_pos_arr  = x_pos_mf->array(mfi);
         auto z_pos_arr  = z_pos_mf->array(mfi);
  
         int k = 0;
         for (int i(il); i <= iu; ++i)
             ofile << x_pos_arr(i,j,k) << ' ' << z_pos_arr(i,j,k) << "\n";
     }
     ofile.close();
 }