ERF_MOSTAverage.H

Go to the documentation of this file.

#ifndef ERF_MOSTAverage_H
#define ERF_MOSTAverage_H
 
#include "AMReX_Gpu.H"
#include "AMReX_FArrayBox.H"
#include "AMReX_MultiFab.H"
#include "AMReX_iMultiFab.H"
#include "AMReX_ParmParse.H"
 
#include "ERF_IndexDefines.H"
#include "ERF_TerrainMetrics.H"
#include "ERF_DataStruct.H"
#include "ERF_EB.H"
#include "ERF_PlanarBoundary.H"
 
class MOSTAverage {
public:
    /**
     * Construct the MOST averaging helper.
     *
     * @param[in] face           orientation of face for MOST averaging
     * @param[in] geom           geometry for all AMR levels
     * @param[in] has_zphys      whether nodal physical-height data are available
     * @param[in] a_pp_prefix    ParmParse prefix
     * @param[in] m_mesh_type    mesh type
     * @param[in] m_terrain_type terrain representation
     * @param[in] zlevels_stag   nominal staggered z levels for all AMR levels
     * @param[in] eb_vec         optional embedded-boundary geometry data
     */
    explicit MOSTAverage (amrex::Orientation face,
                          amrex::Vector<amrex::Geometry>  geom,
                          const bool& has_zphys,
                          std::string a_pp_prefix,
                          const MeshType& m_mesh_type,
                          const TerrainType& m_terrain_type,
                          const amrex::Vector<amrex::Vector<amrex::Real>>& zlevels_stag,
                          const amrex::Vector<const eb_*>& eb_vec = {});
 
    // MOSTAverage() = default;
 
    /**
     * Destroy the MOST averaging helper.
     */
    ~MOSTAverage ()
    {}
 
    /**
     * Default move constructor.
     */
    MOSTAverage (MOSTAverage&&)  noexcept = default;
 
    /**
     * Deleted move-assignment operator.
     *
     * @param[in] other source object
     */
    MOSTAverage& operator=(MOSTAverage&& other)  noexcept = delete;
 
    /**
     * Deleted copy constructor.
     *
     * @param[in] other source object
     */
    MOSTAverage (const MOSTAverage& other) = delete;
 
    /**
     * Deleted copy-assignment operator.
     *
     * @param[in] other source object
     */
    MOSTAverage& operator=(const MOSTAverage& other) = delete;
 
    /**
     * Make MOST-average data structures at one level.
     *
     * @param[in] lev        level index
     * @param[in] vars_old   old-time state and velocity fields
     * @param[in] Theta_prim primitive potential-temperature field
     * @param[in] Qv_prim    primitive water-vapor field
     * @param[in] Qr_prim    primitive rain-water field
     * @param[in] z_phys_nd  nodal physical-height field
     */
    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);
 
    /**
     * Reset pointers to field MultiFabs.
     *
     * @param[in] lev        level index
     * @param[in] vars_old   old-time state and velocity fields
     * @param[in] Theta_prim primitive potential-temperature fields by level
     * @param[in] Qv_prim    primitive water-vapor fields by level
     * @param[in] Qr_prim    primitive rain-water fields by level
     */
    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);
 
    /**
     * Update the rotated fields.
     *
     * @param[in] lev level index
     */
    void set_rotated_fields (const int& lev);
 
    /**
     * Compute number of cells per averaging plane.
     *
     * @param[in] lev level index
     */
    void set_plane_normalization (const int& lev);
 
    /**
     * Compute total embedded-boundary surface area.
     *
     * @param[in] lev level index
     */
    void set_eb_normalization (const int& lev);
 
    /**
     * Compute number of cells in the region average.
     *
     * @param[in] lev level index
     */
    void set_region_normalization (const int& /*lev*/)
    {m_ncell_region = (2 * m_radius + 1) * (2 * m_radius + 1) * (2 * m_radius + 1);}
 
    /**
     * Populate the 2D k-index iMultiFab without terrain.
     *
     * @param[in] lev level index
     */
    void set_k_indices_N (const int& lev);
 
    /**
     * Verify that the lateral face-owned FAB contains the reference stencil.
     *
     * @param[in] lev level index
     *
     * This is a conservative preflight check for the direct regional-average
     * path. It verifies that the reference point plus averaging radius, with
     * allowance for staggered-field neighbors, is available in the owned FAB.
     */
    void validate_lateral_reference_stencil (const int& lev);
 
    /**
     * Populate the 2D k-index iMultiFab with terrain.
     *
     * @param[in] lev level index
     */
    void set_k_indices_T (const int& lev);
 
    /**
     * Populate all 2D normal-index iMultiFabs with terrain.
     *
     * @param[in] lev level index
     */
    void set_norm_indices_T (const int& lev);
 
    /**
     * Populate terrain-aware z positions.
     *
     * @param[in] lev level index
     */
    void set_z_positions_T (const int& lev);
 
    /**
     * Populate embedded-boundary z positions.
     *
     * @param[in] lev level index
     */
    void set_z_positions_EB (const int& lev);
 
    /**
     * Populate terrain-aware normal positions.
     *
     * @param[in] lev level index
     */
    void set_norm_positions_T (const int& lev);
 
    /**
     * Driver for the selected average policy.
     *
     * @param[in] lev level index
     */
    void compute_averages (const int& lev);
 
    /**
     * Fill averages for the plane policy.
     *
     * @param[in] lev level index
     */
    void compute_plane_averages (const int& lev);
 
    /**
     * Fill averages for the point or region policy.
     *
     * @param[in] lev level index
     */
    void compute_region_averages (const int& lev);
 
    /**
     * Fill averages for the embedded-boundary policy.
     *
     * @param[in] lev level index
     */
    void compute_eb_averages (const int& lev);
 
    /**
     * Write k-index data.
     *
     * @param[in] lev level index
     */
    void write_k_indices (const int& lev);
 
    /**
     * Write normal-index data.
     *
     * @param[in] lev level index
     */
    void write_norm_indices (const int& lev);
 
    /**
     * Write XZ planar positions.
     *
     * @param[in] lev level index
     * @param[in] j   y-index of the plane to write
     */
    void write_xz_positions (const int& lev,
                             const int& j);
 
    /**
     * Write averages on the 2D MultiFab.
     *
     * @param[in] lev level index
     */
    void write_averages (const int& lev);
 
    /**
     * Return one 2D average MultiFab.
     *
     * @param[in] lev  level index
     * @param[in] comp average component index
     */
    [[nodiscard]] const amrex::MultiFab* get_average (const int& lev, const int& comp) const { return m_averages[lev][comp].get(); }
 
    /**
     * Return one 2D average MultiFab for modification (used to restore the
     * exponential filter state from a checkpoint).
     *
     * @param[in] lev  level index
     * @param[in] comp average component index
     */
    [[nodiscard]] amrex::MultiFab* get_average (const int& lev, const int& comp) { return m_averages[lev][comp].get(); }
 
    /**
     * Return whether the averages are filtered in time (erf.most.time_average).
     */
    [[nodiscard]] bool do_time_averaging () const { return m_t_avg; }
 
    /**
     * Return the number of average components.
     */
    [[nodiscard]] int get_navg () const { return m_navg; }
 
    /**
     * Return whether the time filter at this level holds meaningful history.
     * Until it does, the next average is taken to be the instantaneous one.
     *
     * @param[in] lev level index
     */
    [[nodiscard]] bool time_avg_is_initialized (const int& lev) const
    {
        return ( m_t_avg && (lev < static_cast<int>(m_t_init.size())) && (m_t_init[lev] == 1) );
    }
 
    /**
     * 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.
     *
     * @param[in] lev level index
     */
    void set_time_avg_initialized (const int& lev)
    {
        if (m_t_avg) { m_t_init[lev] = 1; }
    }
 
    /**
     * 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.
     *
     * @param[in] lev level index
     */
    [[nodiscard]] amrex::Vector<amrex::Real> get_plane_average (const int& lev) const
    {
        if (lev < static_cast<int>(m_plane_average.size())) { return m_plane_average[lev]; }
        return amrex::Vector<amrex::Real>{};
    }
 
    /**
     * 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.
     *
     * @param[in] lev  level index
     * @param[in] pavg filtered plane averages, one per average component
     */
    bool set_plane_average (const int& lev, const amrex::Vector<amrex::Real>& pavg)
    {
        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;
    }
 
    /**
     * Return z_ref, which may be computed from a specified k-index.
     *
     * @param[in] lev level index
     */
    [[nodiscard]] amrex::MultiFab* get_zref (const int& lev) const { return m_zref[lev].get(); }
 
    /**
     * Return the k-index iMultiFab.
     *
     * @param[in] lev level index
     */
    [[nodiscard]] const amrex::iMultiFab* get_k_indices (const int& lev) const { return m_k_indx[lev].get(); }
 
    /**
     * Whether a height lies in the vertical extent [z_lo, z_hi) of a cell.
     *
     * The interval is closed below so that a query height landing exactly on a
     * face (e.g. the default 10 m on a flat mesh with a 10 m first cell) belongs
     * to the cell above it rather than to neither cell.
     *
     * @param[in] z    query height
     * @param[in] z_lo height of the lower face
     * @param[in] z_hi height of the upper face
     */
    AMREX_GPU_HOST_DEVICE AMREX_INLINE
    static bool in_cell_z (const amrex::Real& z,
                           const amrex::Real& z_lo,
                           const amrex::Real& z_hi)
    {
        return ( (z >= z_lo) && (z < z_hi) );
    }
 
    /**
     * Index of the highest cell whose center lies at or below a height above the
     * bottom of the domain, on a column of staggered (face) heights.  This is the
     * stretched-mesh counterpart of floor(zref/dz - 1/2) on a uniform mesh.
     *
     * @param[in] zlevels staggered heights of the column (nz+1 faces)
     * @param[in] zref    height above zlevels[0]
     * @return the cell index, or -1 if zref lies below the first cell center
     */
    static int k_index_below (const amrex::Vector<amrex::Real>& zlevels,
                              const amrex::Real& zref)
    {
        const int nz = static_cast<int>(zlevels.size()) - 1;
        int lk = -1;
        for (int k(0); k < nz; ++k) {
            if (cell_center_height(zlevels, k) <= zref) { lk = k; } else { break; }
        }
        return lk;
    }
 
    /**
     * Height of a cell center above the bottom of the domain, on a column of
     * staggered (face) heights.
     *
     * @param[in] zlevels staggered heights of the column (nz+1 faces)
     * @param[in] k       cell index
     */
    static amrex::Real cell_center_height (const amrex::Vector<amrex::Real>& zlevels,
                                           const int& k)
    {
        return myhalf * (zlevels[k] + zlevels[k+1]) - zlevels[0];
    }
 
    /**
     * Function to compute trilinear interpolation with terrain.
     *
     * @param[in] xp X-position
     * @param[in] yp Y-position
     * @param[in] zp Z-position
     * @param[out] interp_vals Values interpolated
     * @param[in] interp_array Array to interpolate on
     * @param[in] z_arr Physical heights
     * @param[in] plo Problem lower bounds
     * @param[in] dxi Inverse cell size array
     * @param[in] interp_comp Number of components to interpolate
     */
    AMREX_GPU_HOST_DEVICE AMREX_INLINE
    static 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)
    {
        // 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 (in_cell_z(z_target, z_lo, z_hi)) {
                found = true;
                // Weight between the two cell centers that bracket z_target by
                // their physical heights, so the interpolation is linear in z
                // on a stretched or terrain-following column.  A weight taken
                // from the fraction of the containing cell alone (as before)
                // is exact only for equal cell heights: at the 10 m face of a
                // 10 m / 11 m pair it returned the wind at 10.25 m.  A center
                // beyond the nodes searched here is mirrored about the shared
                // face, which gives the equal-height weight at those ends.
                amrex::Real zc = myhalf * (z_lo + z_hi);
                if (z_target >= zc) {
                    amrex::Real zc_up = two * z_hi - zc;
                    if (lk + 2 <= kmax) {
                        amrex::Real z_up = fourth * ( z_arr(i_new,j_new  ,lk+2) + z_arr(i_new+1,j_new  ,lk+2)
                                                  + z_arr(i_new,j_new+1,lk+2) + z_arr(i_new+1,j_new+1,lk+2) );
                        zc_up = myhalf * (z_hi + z_up);
                    }
                    zval = (amrex::Real) (lk + 1) + (z_target - zc) / (zc_up - zc);
                } else {
                    amrex::Real zc_dn = two * z_lo - zc;
                    if (lk >= 1) {
                        amrex::Real z_dn = 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) );
                        zc_dn = myhalf * (z_dn + z_lo);
                    }
                    zval = (amrex::Real) lk + (z_target - zc_dn) / (zc - zc_dn);
                }
                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);
        //       in z, floor(zval) is the upper of the two bracketing cells
        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);
        }
    }
 
    /**
     * Return the number of ghost cells of the 2D data (averages, indices and
     * positions) that can be filled from the field data at this level.
     *
     * With more than one level we must explicitly fill the ghost cells of the 2D
     * data: an isolated fine patch has no neighboring box for a FillBoundary to
     * communicate with, while the fields themselves do carry valid ghost data
     * (filled from the coarse level).  The 2D data are allocated with as many
     * ghost cells as the fields they are built from, but an average reads the
     * fields over m_radius cells in each direction and the tangential velocity
     * magnitude (like the interpolation stencil) reaches one cell beyond that.
     * So a ghost cell of the 2D data can only be computed where the fields carry
     * m_radius+1 ghost cells past it.
     *
     * @param[in] lev level index
     */
    [[nodiscard]] amrex::IntVect get_ng_fill (const int& lev) const;
 
    /**
     * Fill the ghost cells of one planar average, and the valid region of its
     * uncomputed copies when the 3D BoxArray is split in z (see PlanarBoundary).
     *
     * @param[in]     lev level index
     * @param[in,out] mf  planar average to fill
     */
    void fill_planar_boundary (const int& lev, amrex::MultiFab& mf);
 
    /**
     * Fill the ghost cells of one average that cannot be computed from the field
     * data with a zeroth-order extrapolation of the nearest computed value.
     *
     * The averages hold as many ghost cells as the fields, which is more than
     * can be computed from them (see get_ng_fill), so the outermost layers are
     * filled here to leave no average undefined.  This is done before the
     * FillBoundary, which then overwrites whatever is shared with a neighbor.
     *
     * @param[in] lev     level index
     * @param[in] iavg    average component index
     * @param[in] ng_fill number of ghost cells that hold computed data
     */
    void extrap_ghost_cells (const int& lev,
                             const int& iavg,
                             const amrex::IntVect& ng_fill);
 
protected:
 
    // Passed through constructor
    //--------------------------------------------
    amrex::Orientation m_face;                                       // Face to perform averages on
    const amrex::Vector<amrex::Geometry> m_geom;                     // Geometry at each level
    amrex::Vector<amrex::Vector<amrex::MultiFab*>> m_fields;         // Ptr to fields to be averaged
    amrex::Vector<amrex::MultiFab*> m_z_phys_nd;                     // Ptr to terrain height coords
    std::string m_pp_prefix;                                         // ParmParse prefix
    MeshType m_mesh_type;                                            // Mesh type
    TerrainType m_terrain_type;                                      // Terrain type
    amrex::Vector<amrex::Vector<amrex::Real>> m_zlevels_stag;        // Nominal staggered z levels (maxlev)
 
    // General vars for multiple or all policies
    //--------------------------------------------
    int  m_nvar{6};                                                  // 6 fields for U/V/W/T/Qv/Qr
    int  m_navg{9};                                                  // 9 averages for U/V/W/T/Qv/Tv/Umag_XY,Umag_XZ,Umag_YZ
    int  m_maxlev{0};                                                // Total number of levels
    int  m_policy{0};                                                // Policy for type of averaging
    bool m_rotate{false};                                            // Do vector rotations for terrain?
    amrex::Vector<std::unique_ptr<amrex::MultiFab>> m_zref;          // Height above surface for MOST BC
    amrex::Vector<std::unique_ptr<amrex::MultiFab>> m_x_pos;         // Ptr to 2D mf to hold x position (maxlev)
    amrex::Vector<std::unique_ptr<amrex::MultiFab>> m_y_pos;         // Ptr to 2D mf to hold y position (maxlev)
    amrex::Vector<std::unique_ptr<amrex::MultiFab>> m_z_pos;         // Ptr to 2D mf to hold z position (maxlev)
    amrex::Vector<std::unique_ptr<amrex::iMultiFab>> m_i_indx;       // Ptr to 2D imf to hold i indices (maxlev)
    amrex::Vector<std::unique_ptr<amrex::iMultiFab>> m_j_indx;       // Ptr to 2D imf to hold j indices (maxlev)
    amrex::Vector<std::unique_ptr<amrex::iMultiFab>> m_k_indx;       // Ptr to 2D imf to hold k indices (maxlev)
    amrex::Vector<amrex::Vector<std::unique_ptr<amrex::MultiFab>>> m_averages; // Ptr to 2D mf to hold averages (maxlev,navg)
    amrex::Vector<PlanarBoundary> m_planar_bndry;                    // Fills the planar averages on a split BoxArray (maxlev)
    amrex::Vector<amrex::Vector<std::unique_ptr<amrex::MultiFab>>> m_rot_fields; // Rotated field data
 
    // Vars for planar average policy
    //--------------------------------------------
    amrex::Vector<amrex::Vector<int>> m_ncell_plane;                 // Number of cells in plane (maxlev,navg)
    amrex::Vector<amrex::Vector<amrex::Real>> m_plane_average;       // Plane avgs (maxlev,navg)
 
    // Vars for point/region average policy
    //--------------------------------------------
    int m_radius{0};                                                 // Radius around k_index
    int m_ncell_region{1};                                           // Number of cells in local region
    amrex::Vector<int> m_k_in;                                       // Specified k_index for region avg (maxlev)
 
    // Vars for normal vector policy
    //--------------------------------------------
    bool m_interp{false};                                            // Do interpolation on destination?
    bool m_norm_vec{false};                                          // Use normal vector to find IJK?
 
    // Vars for EB averaging policy
    //--------------------------------------------
    amrex::Vector<const eb_*> m_eb_vec;                              // Vector of pointers to EB objects (one per level)
    amrex::Vector<amrex::Vector<amrex::Real>> m_total_bndry_area;   // Total EB surface areas (maxlev, navg)
 
    // Time average w/ exponential filter fun
    //--------------------------------------------
    bool m_t_avg{false};                                             // Flag to do moving average in time
    amrex::Vector<int> m_t_init;                                     // Flag to specify if averages are initialized
    double m_time_window{1.0e-16};                                                // Width of the exp filter function
    amrex::Real m_fact_new, m_fact_old;                              // Time average factors for new and old means
 
    // Surface velocity correction
    //--------------------------------------------
    bool include_subgrid_vel = false;
    amrex::Vector<amrex::Real> m_Vsg;                                // Subgrid velocity scale (Mahrt & Sun 1995 MWR)
 
    // Default values
    //--------------------------------------------
    const amrex::Real zref_default = amrex::Real(10.0);
 
    // Sentinel used to detect whether most.zref was set by the user.  zref_default
    // is itself a legal user value, so "was it specified" cannot be read off the
    // value unless we probe with something unreachable first.  For lateral and
    // z-high faces, most.zref is a wall-relative distance; z-low retains its
    // legacy absolute-coordinate convention. See the parse sites in
    // ERF_MOSTAverage.cpp: they use ParmParse::queryAdd, whose return value is
    // only meaningful on the first parse of a key, and these routines run once
    // per level.
    const amrex::Real zref_sentinel = amrex::Real(-1.0);
};
#endif