MOSTAverage Class Reference

#include <ERF_MOSTAverage.H>

Collaboration diagram for MOSTAverage:

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_average (const int &lev, const int &comp)
bool	do_time_averaging () const
int	get_navg () const
bool	time_avg_is_initialized (const int &lev) const
void	set_time_avg_initialized (const int &lev)
amrex::Vector< amrex::Real >	get_plane_average (const int &lev) const
bool	set_plane_average (const int &lev, const amrex::Vector< amrex::Real > &pavg)
amrex::MultiFab *	get_zref (const int &lev) const
const amrex::iMultiFab *	get_k_indices (const int &lev) const

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

    {}

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:

do_time_averaging()

bool MOSTAverage::do_time_averaging ( ) const

inline

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

{ return m_t_avg; }

Referenced by SurfaceLayer::mac_avg_is_time_averaged().

Here is the caller graph for this function:

get_average() [1/2]

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

inline

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

Parameters

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

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

get_average() [2/2]

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

inline

Return one 2D average MultiFab.

Parameters

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

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

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

Here is the caller graph for this function:

get_k_indices()

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

inline

Return the k-index iMultiFab.

Parameters

[in]

lev

level index

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

get_navg()

int MOSTAverage::get_navg ( ) const

inline

Return the number of average components.

{ return m_navg; }

Referenced by SurfaceLayer::get_num_mac_avg().

Here is the caller graph for this function:

get_plane_average()

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

inline

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

Parameters

[in]

lev

level index

    {
        if (lev < static_cast<int>(m_plane_average.size())) { return m_plane_average[lev]; }
        return amrex::Vector<amrex::Real>{};
    }

Referenced by SurfaceLayer::get_mac_plane_avg().

Here is the caller graph for this function:

get_zref()

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

inline

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

Parameters

[in]

lev

level index

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

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);
 
        // We have just (re)built the average containers at this level, so whatever
        // filter history this level held is gone: m_averages holds bogus_large_value
        // and set_plane_normalization has zeroed m_plane_average.  Note that resize
        // above does not touch existing entries, so this must be set explicitly or a
        // regrid would blend the bogus values into the filtered average.  On restart
        // ReadCheckpointFileSurfaceLayer restores the state and sets this back to 1.
        m_t_init[lev] = 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);
    // See zref_sentinel in ERF_MOSTAverage.H: probe with the sentinel and test the
    // value, since queryAdd's return value stops meaning "user specified" after the
    // first parse of the key and this routine runs once per level.
    Real zref_tmp = zref_sentinel;
    pp.queryAdd("most.zref",zref_tmp);
    bool read_z = (zref_tmp > Real(0));
    if (!read_z) { zref_tmp = zref_default; }
    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);
    // See zref_sentinel in ERF_MOSTAverage.H: probe with the sentinel and test the
    // value, since queryAdd's return value stops meaning "user specified" after the
    // first parse of the key and this routine runs once per level.
    Real zref_tmp = zref_sentinel;
    pp.queryAdd("most.zref",zref_tmp);
    bool read_z = (zref_tmp > Real(0));
    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) {
        zref_tmp = zref_default;
        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);
    // See zref_sentinel in ERF_MOSTAverage.H: probe with the sentinel and test the
    // value, since queryAdd's return value stops meaning "user specified" after the
    // first parse of the key and this routine runs once per level.
    Real zref_tmp = zref_sentinel;
    pp.queryAdd("most.zref",zref_tmp);
    bool read_zref = (zref_tmp > Real(0));
    if (!read_zref) {
        zref_tmp = zref_default;
        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);
    // See zref_sentinel in ERF_MOSTAverage.H: probe with the sentinel and test the
    // value, since queryAdd's return value stops meaning "user specified" after the
    // first parse of the key and this routine runs once per level.
    Real zref_tmp = zref_sentinel;
    pp.queryAdd("most.zref",zref_tmp);
    bool read_zref = (zref_tmp > Real(0));
    if (!read_zref) {
        zref_tmp = zref_default;
        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_average()

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

inline

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

Parameters

[in]	lev	level index
[in]	pavg	filtered plane averages, one per average component

    {
        const int n_have = (lev < static_cast<int>(m_plane_average.size())) ?
                            static_cast<int>(m_plane_average[lev].size()) : 0;
        if (n_have != static_cast<int>(pavg.size())) { return false; }
        if (n_have > 0) { m_plane_average[lev] = pavg; }
        return true;
    }

Referenced by SurfaceLayer::set_mac_plane_avg().

Here is the caller graph for this function:

set_plane_normalization()

void MOSTAverage::set_plane_normalization

(

const int &

lev

)

Compute number of cells per averaging plane.

Parameters

[in]

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

    {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_time_avg_initialized()

void MOSTAverage::set_time_avg_initialized ( const int &  lev )

inline

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

Parameters

[in]

lev

level index

    {
        if (m_t_avg) { m_t_init[lev] = 1; }
    }

Referenced by SurfaceLayer::set_mac_avg_initialized().

Here is the caller graph for this function:

set_z_positions_EB()

void MOSTAverage::set_z_positions_EB

(

const int &

lev

)

Populate embedded-boundary z positions.

Parameters

[in]

lev

level index

Function to set K indices for EB.

Parameters

[in]

lev

Current level

{
    // See zref_sentinel in ERF_MOSTAverage.H: probe with the sentinel and test the
    // value, since queryAdd's return value stops meaning "user specified" after the
    // first parse of the key and this routine runs once per level.
    Real zref_tmp = zref_sentinel;
    ParmParse pp(m_pp_prefix);
    pp.queryAdd("most.zref",zref_tmp);
    bool read_z = (zref_tmp > Real(0));
    if (!read_z) { zref_tmp = zref_default; }
 
    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);
    // See zref_sentinel in ERF_MOSTAverage.H: probe with the sentinel and test the
    // value, since queryAdd's return value stops meaning "user specified" after the
    // first parse of the key and this routine runs once per level.
    Real zref_tmp = zref_sentinel;
    pp.queryAdd("most.zref",zref_tmp);
    bool read_zref = (zref_tmp > Real(0));
    if (!read_zref) {
        zref_tmp = zref_default;
        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:

time_avg_is_initialized()

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

inline

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

Parameters

[in]

lev

level index

    {
        return ( m_t_avg && (lev < static_cast<int>(m_t_init.size())) && (m_t_init[lev] == 1) );
    }

Referenced by SurfaceLayer::mac_avg_is_initialized().

Here is the caller graph for this function:

trilinear_interp_T()

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

inlinestatic

Function to compute trilinear interpolation with terrain.

Parameters

[in]	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();
}

Here is the call graph for this function:

Member Data Documentation

include_subgrid_vel

bool MOSTAverage::include_subgrid_vel = false

protected

Referenced by make_MOSTAverage_at_level().

m_averages

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

protected

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

m_eb_vec

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

protected

Referenced by compute_eb_averages(), and set_eb_normalization().

m_fact_new

amrex::Real MOSTAverage::m_fact_new

protected

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

m_fact_old

amrex::Real MOSTAverage::m_fact_old

protected

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

m_fields

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

protected

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

m_geom

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

protected

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

m_i_indx

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

protected

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

m_interp

bool MOSTAverage::m_interp {false}

protected

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

m_j_indx

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

protected

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

m_k_in

amrex::Vector<int> MOSTAverage::m_k_in

protected

Referenced by set_k_indices_N(), and set_k_indices_T().

m_k_indx

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

protected

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

m_maxlev

int MOSTAverage::m_maxlev {0}

protected

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

m_mesh_type

MeshType MOSTAverage::m_mesh_type

protected

m_navg

int MOSTAverage::m_navg {6}

protected

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

m_ncell_plane

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

protected

Referenced by compute_plane_averages(), and set_plane_normalization().

m_ncell_region

int MOSTAverage::m_ncell_region {1}

protected

Referenced by compute_region_averages(), and set_region_normalization().

m_norm_vec

bool MOSTAverage::m_norm_vec {false}

protected

Referenced by make_MOSTAverage_at_level().

m_nvar

int MOSTAverage::m_nvar {6}

protected

Referenced by make_MOSTAverage_at_level().

m_plane_average

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

protected

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

m_policy

int MOSTAverage::m_policy {0}

protected

Referenced by compute_averages(), and make_MOSTAverage_at_level().

m_pp_prefix

std::string MOSTAverage::m_pp_prefix

protected

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

m_radius

int MOSTAverage::m_radius {0}

protected

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

m_rot_fields

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

protected

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

m_rotate

bool MOSTAverage::m_rotate {false}

protected

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

m_t_avg

bool MOSTAverage::m_t_avg {false}

protected

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

m_t_init

amrex::Vector<int> MOSTAverage::m_t_init

protected

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

m_terrain_type

TerrainType MOSTAverage::m_terrain_type

protected

Referenced by make_MOSTAverage_at_level().

m_time_window

double MOSTAverage::m_time_window {1.0e-16}

protected

Referenced by make_MOSTAverage_at_level().

m_total_bndry_area

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

protected

Referenced by compute_eb_averages(), and set_eb_normalization().

m_Vsg

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

protected

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

m_x_pos

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

protected

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

m_y_pos

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

protected

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

m_z_phys_nd

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

protected

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

m_z_pos

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

protected

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

m_zref

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

protected

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

zref_default

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

protected

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

zref_sentinel

const amrex::Real MOSTAverage::zref_sentinel = amrex::Real(-1.0)

protected

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

---

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

• Source/BoundaryConditions/ERF_MOSTAverage.H
• Source/BoundaryConditions/ERF_MOSTAverage.cpp