ERF_TurbPertStruct.H

Go to the documentation of this file.

#ifndef ERF_TURB_PERT_STRUCT_H_
#define ERF_TURB_PERT_STRUCT_H_
 
#include <ERF_DataStruct.H>
#include <AMReX_MultiFabUtil.H>
#include <ERF_TileNoZ.H>
#include <time.h>
/**
 * Container holding quantities related to turbulent perturbation parameters
 */
 
/*  The general rule of thumb is to create a perturbation box size of 1/8th of the boundary layer height.
    The boundary layer height can't be the height of the domain. The length and width of the box should
    be twice the height of the box. If meandering flow is present, the width of the box should be take
    the angle of the inflow into consideration.
*/
 
AMREX_ENUM(PerturbationType,
    Source, Direct, CPM, CPM_W, None
);
 
/**
 * @brief Runtime state and operations for turbulent perturbation forcing.
 */
struct TurbulentPerturbation {
 
  public:
 
    /**
     * @brief Destroy turbulent perturbation state.
     */
    ~TurbulentPerturbation () {}
 
    /**
     * @brief Store the integer perturbation-type selector for one AMR level.
     * @param lev AMR level index.
     * @param pert_type Perturbation type selected for this level.
     * @param max_level Maximum AMR level configured for the run.
     */
    void init_tpi_type (const int lev,
                        const PerturbationType& pert_type,
                        const int max_level)
    {
        if (pt_type.size() < max_level + 1) {
            pt_type.resize(max_level + 1, -1);
        }
 
        if (pert_type == PerturbationType::Source) {
            pt_type[lev] = 0;
        } else if (pert_type == PerturbationType::Direct) {
            pt_type[lev] = 1;
        } else if (pert_type == PerturbationType::CPM) {
            pt_type[lev] = 2;
        } else if (pert_type == PerturbationType::CPM_W) {
            pt_type[lev] = 3;
        } else {
            pt_type[lev] = -1;
        }
    }
 
    /**
     * @brief Initialize turbulent perturbation regions and per-box storage.
     * @param lev AMR level index.
     * @param subdomains_lev Rectangular level subdomains used to define perturbation regions.
     * @param dx Cell spacing for this level.
     * @param ba State BoxArray used for perturbation-cell storage.
     * @param dm DistributionMapping for perturbation-cell storage.
     * @param ngrow_state Number of ghost cells for perturbation-cell storage.
     * @param pp_prefix ParmParse prefix for the ERF input namespace.
     * @param refRatio Refinement ratios between AMR levels.
     * @param max_level Maximum AMR level configured for the run.
     */
    void init_tpi (const int lev,
                   const amrex::Vector<amrex::BoxArray>& subdomains_lev,
                   const amrex::GpuArray<amrex::Real,3> dx,
                   const amrex::BoxArray& ba,
                   const amrex::DistributionMapping& dm,
                   const int ngrow_state,
                   std::string pp_prefix,
                   const amrex::Vector<amrex::IntVect> refRatio,
                   const int max_level)
 
    {
        ref_ratio = refRatio;
 
        amrex::ParmParse pp(pp_prefix);
 
        // Reading inputs, and placing assertion for the perturbation inflow to work
        pp.getarr("perturbation_box_dims",tpi_boxDim);
        pp.getarr("perturbation_direction",tpi_direction);
        pp.get("perturbation_layers",tpi_layers);
        pp.get("perturbation_offset",tpi_offset);
 
        tpi_nonDim = zero;
        pp.query("perturbation_nondimensional",tpi_nonDim);
 
        tpi_Tinf = amrex::Real(300.);
        pp.query("perturbation_T_infinity",tpi_Tinf);
 
        tpi_Ti = zero;
        pp.query("perturbation_T_intensity",tpi_Ti);
 
        input_Ug = zero;
        pp.query("perturbation_Ug",input_Ug);
 
        input_w_amp = zero;
        pp.query("perturbation_w_amp",input_w_amp);
 
        perturbation_klo = 0;
        bool have_klo = pp.query("perturbation_klo", perturbation_klo);
 
        perturbation_khi = 0;
        bool have_khi = pp.query("perturbation_khi", perturbation_khi);
 
        // Check variables message
        if (tpi_offset < 0) { amrex::Abort("Please provide a valid inflow cell offset value for perturbation region (ie. 0-5)"); }
        if (tpi_layers < 0) { amrex::Abort("Please provide a valid perturbation layer value (ie. 3-5)"); }
        if (tpi_nonDim < zero) { amrex::Abort("Please provide a valid nondimensional number (ie. Ri = amrex::Real(0.042))"); }
        for (int i = 0; i < tpi_boxDim.size(); i++) {
            if (tpi_boxDim[i] < 3) { amrex::Abort("Please provide valid dimensions for perturbation boxes."); }
        }
        if (input_Ug   < zero) { amrex::Abort("Please provide a valid geostrophic wind speed (ie. Ug = amrex::Real(10.0) m/s)"); }
        if (tpi_Tinf   < zero) { amrex::Abort("Please provide a valid ambient temperature value (ie. T_0 = T_infty)"); }
        if (tpi_Ti     < zero) { amrex::Abort("Please provide a valid temperature intensity value (ie. 0-one)"); }
 
        // Create a temporary box list to accumulate all the perturbation regions after box modification
        amrex::BoxList tmp_bl;
 
        // boxSize for individual boxes
        amrex::IntVect boxSize(tpi_boxDim[0],tpi_boxDim[1],tpi_boxDim[2]);
 
        if (tpi_direction[2] || tpi_direction[5]) { amrex::Abort("Currently not supporting z-direction flow perturbation"); }
 
        for (int isub = 0; isub < subdomains_lev.size(); ++isub) {
            const amrex::BoxArray& subdomain = subdomains_lev[isub];
            amrex::Box subdomain_box(subdomain.minimalBox());
 
            if (subdomain_box.numPts() != subdomain.numPts()) {
                amrex::Abort("Turbulent perturbations require rectangular subdomains. "
                             "Level " + std::to_string(lev) +
                             ", subdomain " + std::to_string(isub) +
                             " is not a rectangular region fully covered by grids.");
            }
 
            const amrex::IntVect& valid_box_lo = subdomain_box.smallEnd();
            const amrex::IntVect& valid_box_hi = subdomain_box.bigEnd();
 
            // default perturbation region to be the entire z extent if the user does not specify bounds
            if (!have_klo) { perturbation_klo = valid_box_lo[2]; }
            if (!have_khi) { perturbation_khi = valid_box_hi[2]; }
 
            // Creating perturbation regions and initializing with generic size.
            amrex::Box lo_x_bx(amrex::IntVect(0), amrex::IntVect(1), amrex::IntVect(0));
            amrex::Box hi_x_bx(amrex::IntVect(0), amrex::IntVect(1), amrex::IntVect(0));
            amrex::Box lo_y_bx(amrex::IntVect(0), amrex::IntVect(1), amrex::IntVect(0));
            amrex::Box hi_y_bx(amrex::IntVect(0), amrex::IntVect(1), amrex::IntVect(0));
 
            // Starting logic to set the size of the perturbation region(s)
            //amrex::PrintToFile("BoxPerturbationOutput") << "Setting perturbation region in:";
            // ***** X-direction perturbation *****
            if (tpi_direction[0]) { // West
                lo_x_bx.setSmall(amrex::IntVect(valid_box_lo[0]+tpi_offset, valid_box_lo[1]+tpi_direction[1]*tpi_offset, perturbation_klo));
                lo_x_bx.setBig  (amrex::IntVect(valid_box_lo[0]+(tpi_layers*tpi_boxDim[0]-1)+tpi_offset, valid_box_hi[1]-(tpi_direction[4]*tpi_offset), perturbation_khi));
                amrex::PrintToFile("BoxPerturbationOutput") << " West face";
            }
 
            if (tpi_direction[3]) { // East
                hi_x_bx.setSmall(amrex::IntVect(valid_box_hi[0]-((tpi_layers*tpi_boxDim[0]-1)+tpi_offset), valid_box_lo[1]+tpi_direction[1]*tpi_offset, perturbation_klo));
                hi_x_bx.setBig(amrex::IntVect(valid_box_hi[0]-tpi_offset, valid_box_hi[1]-(tpi_direction[4]*tpi_offset), perturbation_khi));
                amrex::PrintToFile("BoxPerturbationOutput") << " East face";
            }
 
            // ***** Y-direction Perturbation *****
            if (tpi_direction[1]) { // North
                lo_y_bx.setSmall(amrex::IntVect(valid_box_lo[0]+tpi_direction[0]*tpi_offset, valid_box_lo[1]+tpi_offset, perturbation_klo));
                lo_y_bx.setBig(amrex::IntVect(valid_box_hi[0]-tpi_direction[3]*tpi_offset, valid_box_lo[1]+((tpi_layers*tpi_boxDim[1])-1)+tpi_offset, perturbation_khi));
                amrex::PrintToFile("BoxPerturbationOutput") << " North face";
            }
 
            if (tpi_direction[4]) { // South
                hi_y_bx.setSmall(amrex::IntVect(valid_box_lo[0]+tpi_direction[0]*tpi_offset, valid_box_hi[1]-((tpi_layers*tpi_boxDim[1]-1)+tpi_offset), perturbation_klo));
                hi_y_bx.setBig(amrex::IntVect(valid_box_hi[0]-tpi_direction[3]*tpi_offset, valid_box_hi[1]-tpi_offset, perturbation_khi));
                amrex::PrintToFile("BoxPerturbationOutput") << " South face";
            }
 
            // Performing box union for intersecting perturbation regions to avoid overlapping sections (double counting at corners)
            if (tpi_direction[0] && tpi_direction[1]) {  // Reshaping South smallEnd
                amrex::Box lo_x_lo_y_u = lo_x_bx & lo_y_bx;
                lo_y_bx.setSmall(amrex::IntVect(lo_x_lo_y_u.bigEnd(0)+1, lo_x_lo_y_u.smallEnd(1), lo_x_lo_y_u.smallEnd(2)));
            }
 
            if (tpi_direction[3] && tpi_direction[1]) { // Reshaping South bigEnd
                 amrex::Box hi_x_lo_y_u = hi_x_bx & lo_y_bx;
                 lo_y_bx.setBig(amrex::IntVect(hi_x_lo_y_u.smallEnd(0)-1, hi_x_lo_y_u.bigEnd(1), hi_x_lo_y_u.bigEnd(2)));
            }
 
            if (tpi_direction[0] && tpi_direction[4]) {  // Reshaping North smallEnd
                amrex::Box lo_x_hi_y_u = lo_x_bx & hi_y_bx;
                hi_y_bx.setSmall(amrex::IntVect(lo_x_hi_y_u.bigEnd(0)+1, lo_x_hi_y_u.smallEnd(1), lo_x_hi_y_u.smallEnd(2)));
            }
 
            if (tpi_direction[3] && tpi_direction[4]) { // Reshaping North bigEnd
                 amrex::Box hi_x_hi_y_u = hi_x_bx & hi_y_bx;
                 hi_y_bx.setBig(amrex::IntVect(hi_x_hi_y_u.smallEnd(0)-1, hi_x_hi_y_u.bigEnd(1), hi_x_hi_y_u.bigEnd(2)));
            }
 
            // Creating structure box array for conserved quantity
            if (tpi_direction[0]) { tmp_bl.push_back(lo_x_bx); }
            if (tpi_direction[1]) { tmp_bl.push_back(lo_y_bx); }
            if (tpi_direction[3]) { tmp_bl.push_back(hi_x_bx); }
            if (tpi_direction[4]) { tmp_bl.push_back(hi_y_bx); }
        }
 
        //amrex::PrintToFile("BoxPerturbationOutput") << "\nBoxList: " << tmp_bl << "\n";
        amrex::BoxArray tmp_ba(tmp_bl);
        tmp_ba.maxSize(boxSize);
 
        const int num_levels = max_level + 1;
        if (pb_ba.size() < num_levels) {
            pb_ba.resize(num_levels);
            pb_mag.resize(num_levels);
            pb_dir.resize(num_levels);
            pb_netZero.resize(num_levels);
            pb_interval.resize(num_levels);
            pb_local_etime.resize(num_levels);
            pb_amp.resize(num_levels);
            pb_cell.resize(num_levels);
            tpi_Lpb.resize(num_levels);
            tpi_Wpb.resize(num_levels);
            tpi_Hpb.resize(num_levels);
            tpi_lref.resize(num_levels);
        }
 
        pb_ba[lev] = tmp_ba;
 
        // Initializing mean magnitude and direction vectors
        pb_mag[lev].resize(pb_ba[lev].size(), zero);
        pb_dir[lev].resize(pb_ba[lev].size(), zero);
        pb_netZero[lev].resize(pb_ba[lev].size(), zero);
 
        // Set size of vector and initialize
        pb_interval[lev].resize(pb_ba[lev].size(), -one);
        pb_local_etime[lev].resize(pb_ba[lev].size(), zero);
        pb_amp[lev].resize(pb_ba[lev].size(), zero);
 
        // Creating data array for perturbation amplitude storage
        if (pt_type[lev] == 3) { // CPM_W converts to the k-face ba.
            pb_cell[lev].define(convert(ba, amrex::IntVect(0,0,1)), dm, 1, ngrow_state);
        } else {
            pb_cell[lev].define(ba, dm, 1, ngrow_state); // this is the only place ba is used. Maybe we can print here to determine what's valid...
        }
        pb_cell[lev].setVal(0.);
 
        // Computing perturbation reference length
        tpi_Lpb[lev] = tpi_boxDim[0]*dx[0];
        tpi_Wpb[lev] = tpi_boxDim[1]*dx[1];
        tpi_Hpb[lev] = tpi_boxDim[2]*dx[2];
        tpi_lref[lev] = std::sqrt(tpi_Lpb[lev]*tpi_Lpb[lev] + tpi_Wpb[lev]*tpi_Wpb[lev]);
 
        tpi_pert_adjust = zero;
        tpi_net_buoyant = zero;
 
        /*
        // Function check point message
        amrex::PrintToFile("BoxPerturbationOutput") << "perturbation_box_dims: "
                                                    << tpi_boxDim[0] << " "
                                                    << tpi_boxDim[1] << " "
                                                    << tpi_boxDim[2] << "\n";
        amrex::PrintToFile("BoxPerturbationOutput") << "perturbation_direction: "
                                                    << tpi_direction[0] << " "
                                                    << tpi_direction[1] << " "
                                                    << tpi_direction[2] << "\n\n";
        amrex::PrintToFile("BoxPerturbationOutput") << "perturbation_layers: " << tpi_layers << "\n";
        amrex::PrintToFile("BoxPerturbationOutput") << "perturbation_offset: " << tpi_offset << "\n\n";
        amrex::PrintToFile("BoxPerturbationOutput") << "perturbation_nondimensional: " << tpi_nonDim << "\n";
        amrex::PrintToFile("BoxPerturbationOutput") << "perturbation_T_infinity: " << tpi_Tinf << "\n";
        amrex::PrintToFile("BoxPerturbationOutput") << "perturbation_T_intensity: " << tpi_Ti << "\n";
        amrex::PrintToFile("BoxPerturbationOutput") << "Reference length per box = " << tpi_lref[lev] << "\n\n";
        amrex::PrintToFile("BoxPerturbationOutput") << "Turbulent perturbation BoxArray:\n" << pb_ba[lev] << "\n";
        */
    }
 
 
    /**
     * @brief Update perturbation amplitudes and intervals when each box is ready.
     * @param lev AMR level index.
     * @param dt Time step for elapsed-time accumulation.
     * @param mf_xvel X-velocity field used to compute box mean velocities.
     * @param mf_yvel Y-velocity field used to compute box mean velocities.
     * @param mf_cons Conserved state used for iteration layout.
     */
    void calc_tpi_update (const int lev,
                          const double dt,
                          amrex::MultiFab& mf_xvel,
                          amrex::MultiFab& mf_yvel,
                          amrex::MultiFab& mf_cons)
    {
        // Resetting the net buoyant force value
        tpi_net_buoyant = zero;
 
        // Setting random number generator for update interval
        srand( (unsigned) time(NULL) );
 
        auto m_ixtype = mf_cons.boxArray().ixType(); // safety step
        if (pt_type[lev] == 3) { m_ixtype = amrex::IndexType(amrex::IntVect(0,0,1)); }
 
        // Seed the random generator at 1024UL for regression testing
        int fix_random_seed = 0;
        amrex::ParmParse pp("erf");
        pp.query("fix_random_seed", fix_random_seed);
        if (fix_random_seed) {
            // We need this one for the ParalleForRNG used in calc_tpi
            amrex::InitRandom(1024UL, amrex::ParallelDescriptor::NProcs(), 1024UL);
 
            // We need this one for the RandomReal below
            srand(1024UL);
        }
 
        for (int boxIdx = 0; boxIdx < pb_ba[lev].size(); boxIdx++) {
 
            bool update_box = true; // initialize flag
            // Check if the local elapsed time is greater than the update interval and don't go into boxes the rank doesn't own
            if (pt_type[lev] == 2 || pt_type[lev] == 3) {
                // for CPM, perturbation boxes/cells refresh after boxes have advected box width/length * num_layers (advective time scale)
                update_box = ((pb_local_etime[lev][boxIdx] >= pb_interval[lev][boxIdx]*tpi_layers) && pb_local_etime[lev][boxIdx] != -one);
            } else {
                update_box = ( pb_local_etime[lev][boxIdx] >= pb_interval[lev][boxIdx] && pb_local_etime[lev][boxIdx] != -one );
            }
            if ( update_box ) {
 
                // Compute mean velocity of each perturbation box
                calc_tpi_meanMag_perBox(boxIdx, lev, mf_cons, mf_xvel, mf_yvel);
 
                // Only the rank owning the box will be able to access the storage location
                // Done for parallelism to avoid Inf being stored in array
                if (pb_mag[lev][boxIdx] !=zero) {
                    double interval = 0.0;
                    if (pt_type[lev] == 2 || pt_type[lev] == 3) {
                        //  Wind direction correction for angled wind
                        amrex::Real wind_direction = pb_dir[lev][boxIdx];
                        if (wind_direction > PI / 4) { wind_direction = PI / 2 - wind_direction; }
                        // CPM only cares about the side length of a box, min call maintains flexibility.
                        interval = one / std::cos(wind_direction) * std::min(tpi_Lpb[lev], tpi_Wpb[lev]) / pb_mag[lev][boxIdx];
                    } else {
                        interval = tpi_lref[lev] / pb_mag[lev][boxIdx];
                    }
                    pb_interval[lev][boxIdx] = RandomReal(static_cast<amrex::Real>(0.9*interval),
                                                          static_cast<amrex::Real>(1.1*interval)); // 10% variation
 
                    // Reset local elapsed time
                    pb_local_etime[lev][boxIdx] = zero;
                } else {
                    // this box is not on this rank, we shouldn't enter it again.
                    // Technically, all boxes are looped through on the very first step of a simulation and this is when this is set
                    pb_local_etime[lev][boxIdx] = -one;
                }
 
                // Trigger amplitude calculation per perturbation box
                calc_tpi_amp(lev, boxIdx, pb_interval[lev][boxIdx]);
 
                // Trigger random amplitude storage per cell within perturbation box
                pseudoRandomPert(boxIdx, lev, m_ixtype);
 
            } else {
                // set perturbation amplitudes to 0 for CPM. A little inefficient but leverages as much as existing code as possible.
                if (pt_type[lev] == 2 || pt_type[lev] == 3) { zero_amp(boxIdx, lev, m_ixtype); }
 
                // Increase by timestep of level 0 (but only if the box is owned by this rank)
                if (pb_local_etime[lev][boxIdx] != -one) { pb_local_etime[lev][boxIdx] += dt; }
            } // if
 
            if (pt_type[lev] < 2) { // box perturbation method only
                // Per iteration operation of net-zero buoyant force check
                if (pb_mag[lev][boxIdx] !=zero) netZeroBuoyantAdd(boxIdx, lev);
                tpi_net_buoyant += pb_netZero[lev][boxIdx];
            }
        } // for
 
        if (pt_type[lev] < 2) { // box perturbation method only
            // Normalizing the adjustment based on how many boxes there are
            // the values within the array is already normalized by the number
            // of cells within each box
            tpi_pert_adjust = tpi_net_buoyant / (amrex::Real) pb_ba[lev].size();
 
            // Per iteration operation of net-zero buoyant force adjustment
            for (int boxIdx = 0; boxIdx < pb_ba[lev].size(); boxIdx++) {
                if (pb_mag[lev][boxIdx] !=zero) netZeroBuoyantAdjust(boxIdx, lev);
            }
        }
    }
 
    /**
     * @brief Apply stored turbulent perturbations to a source or state array.
     * @param lev AMR level index.
     * @param vbx Valid box over which perturbations may be applied.
     * @param comp Component to modify, or -1 for vertical velocity perturbations.
     * @param m_ixtype Index type of the destination array.
     * @param src_arr Destination array receiving perturbation increments.
     * @param pert_cell Stored perturbation amplitudes.
     */
    void apply_tpi (const int& lev,
                    const amrex::Box& vbx,                      // box union from upper level
                    const int& comp,                            // Component to modify
                    const amrex::IndexType& m_ixtype,           // IntVect type of src_arr
                    const amrex::Array4<amrex::Real>& src_arr,  // Array to apply perturbation
                    const amrex::Array4<amrex::Real const>& pert_cell)
    {
        for (int boxIdx = 0; boxIdx < pb_ba[lev].size(); boxIdx++) {
            amrex::Box pbx = amrex::convert(pb_ba[lev][boxIdx], m_ixtype);
            if (pt_type[lev] == 3) { pbx.setBig(2, pbx.bigEnd(2) - 1); } // prevent double counting after converting to k-faces
            amrex::Box ubx = pbx & vbx;
            if (ubx.ok()) {
                if (comp == -1) { // vertical velocity perturbations
                    ParallelFor(ubx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
                        src_arr(i,j,k) += pert_cell(i,j,k);
                    });
                } else {
                    ParallelFor(ubx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
                        src_arr(i,j,k,comp) += pert_cell(i,j,k);
 
                        // For box region debug only
                        #ifdef INDEX_PERTURB
                        src_arr(i,j,k,comp) = (amrex::Real) (boxIdx + 5.);
                        #endif
                    });
                }
            }
        }
    }
 
    /**
     * @brief Compute the perturbation amplitude for one perturbation box.
     * @param lev AMR level index.
     * @param boxIdx Perturbation-box index on the level.
     * @param interval Current perturbation update interval [s].
     */
    void calc_tpi_amp (const int& lev,
                       const int& boxIdx,
                       const double& interval)
    {
        pb_amp[lev][boxIdx] = zero; // Safety step
        if (pt_type[lev] == 2) { // CPM
            amrex::Real cp = 1004; // specific heat of air [J/(kg K)]
            amrex::Real Ec = amrex::Real(0.2); // Eckert number
            pb_amp[lev][boxIdx] = (input_Ug * input_Ug) / (Ec * cp);
        } else if (pt_type[lev] == 3) {
            pb_amp[lev][boxIdx] = input_w_amp;
        } else { // box perturbation
            amrex::Real Um = pb_mag[lev][boxIdx];
            amrex::Real beta = one/tpi_Tinf; // Thermal expansion coefficient
 
            // Pseudo Random temperature the ignores scale when mechanically tripping turbulence
            amrex::Real g = CONST_GRAV;
            // get total refinement ratio on each level
            int total_ref_ratio = 1;
            for (int level = lev; level >= 1; level--) {
                total_ref_ratio *= ref_ratio[level-1][2];
            }
            // calculation needs to be scale-aware since the formulation relies on the physical size of the box
            if (tpi_Ti > zero) g = (tpi_nonDim * Um * Um) / (tpi_Ti * tpi_Hpb[lev]) * 1 / total_ref_ratio;
 
            // Ma and Senocak (2023) Eq. 8, solving for delta phi
            pb_amp[lev][boxIdx] = (tpi_nonDim * Um * Um) / (g * beta * tpi_Hpb[lev]) * 1 / total_ref_ratio;
 
            if (pt_type[lev] == 0) {
                // Performing this step converts the perturbation proportionality into
                // the forcing term
                // Ma & Senocak (2023) Eq. 7
                pb_amp[lev][boxIdx] /= static_cast<amrex::Real>(interval);
            }
        }
    }
 
    /**
     * @brief Assign pseudo-random perturbations to cells in one perturbation box.
     * @param boxIdx Perturbation-box index on the level.
     * @param lev AMR level index.
     * @param m_ixtype Index type used to align perturbation storage.
     */
    void pseudoRandomPert (const int& boxIdx,
                           const int& lev,
                           const amrex::IndexType& m_ixtype)
    {
        for (amrex::MFIter mfi(pb_cell[lev],TileNoZ()); mfi.isValid(); ++mfi) {
            amrex::Box vbx = mfi.validbox();
            amrex::Box pbx = amrex::convert(pb_ba[lev][boxIdx], m_ixtype);
            if (pt_type[lev] == 3) { pbx.setBig(2, pbx.bigEnd(2) - 1); } // prevent double counting after converting to k-faces
            amrex::Box ubx = pbx & vbx;
            if (ubx.ok()) {
                amrex::Real amp_copy = pb_amp[lev][boxIdx];
                const amrex::Array4<amrex::Real>& pert_cell = pb_cell[lev].array(mfi);
 
                if (pt_type[lev] == 2  || pt_type[lev] == 3) { // CPM
                    amrex::Real rand_number_const = RandomReal(-one, one);
                    ParallelFor(ubx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
                        pert_cell(i,j,k) = rand_number_const * amp_copy;
                    });
                } else {
                    ParallelForRNG(ubx, [=] AMREX_GPU_DEVICE(int i, int j, int k, const amrex::RandomEngine& engine) noexcept {
                        amrex::Real rand_double = amrex::Random(engine);
                        pert_cell(i,j,k) = (rand_double*two - one) * amp_copy;
                    });
                }
            }
        }
    }
 
    /**
     * @brief Reset CPM perturbation amplitudes in one perturbation box to zero.
     * @param boxIdx Perturbation-box index on the level.
     * @param lev AMR level index.
     * @param m_ixtype Index type used to align perturbation storage.
     */
    void zero_amp (const int& boxIdx,
                   const int& lev,
                   const amrex::IndexType& m_ixtype)
    {
 
        for (amrex::MFIter mfi(pb_cell[lev],TileNoZ()); mfi.isValid(); ++mfi) {
            amrex::Box vbx = mfi.validbox();
            amrex::Box pbx = amrex::convert(pb_ba[lev][boxIdx], m_ixtype);
            if (pt_type[lev] == 3) { pbx.setBig(2, pbx.bigEnd(2) - 1); } // prevent double counting after converting to k-faces
            amrex::Box ubx = pbx & vbx;
            if (ubx.ok()) {
                const amrex::Array4<amrex::Real>& pert_cell = pb_cell[lev].array(mfi);
                ParallelFor(ubx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
                    pert_cell(i,j,k) = zero;
                });
            }
        }
    }
 
    /**
     * @brief Accumulate the mean perturbation contribution for net-zero buoyancy enforcement.
     * @param boxIdx Perturbation-box index on the level.
     * @param lev AMR level index.
     */
    void netZeroBuoyantAdd (const int& boxIdx,
                            const int& lev)
    {
        // Creating local copy of PB box array and magnitude
        const amrex::BoxArray m_pb_ba = pb_ba[lev];
        amrex::Real* m_pb_netZero = pb_netZero[lev].data();
 
        // Create device array for summation
        amrex::Vector<amrex::Real> avg_h(1,zero);
        amrex::Gpu::DeviceVector<amrex::Real> avg_d(1,zero);
        amrex::Real* avg = avg_d.data();
 
        // Iterates through the cells of each box and sum the white noise perturbation
        for (amrex::MFIter mfi(pb_cell[lev], TileNoZ()) ; mfi.isValid(); ++mfi) {
            const amrex::Box& vbx = mfi.validbox();
            amrex::Box pbx = amrex::convert(m_pb_ba[boxIdx], vbx.ixType());
            amrex::Box ubx = pbx & vbx;
            if (ubx.ok()) {
                const amrex::Array4<const amrex::Real>& pert_cell = pb_cell[lev].const_array(mfi);
                amrex::Real norm = one / static_cast<amrex::Real>(ubx.numPts());
                ParallelFor(amrex::Gpu::KernelInfo().setReduction(true), ubx, [=]
                AMREX_GPU_DEVICE(int i, int j, int k, amrex::Gpu::Handler const& handler) noexcept {
                    amrex::Gpu::deviceReduceSum(&avg[0], pert_cell(i,j,k)*norm, handler);
                });
                amrex::Gpu::copy(amrex::Gpu::deviceToHost, avg_d.begin(), avg_d.end(), avg_h.begin());
 
                // Assigning onto storage array
                m_pb_netZero[boxIdx] = avg_h[0];
            }
        }
    }
 
    /**
     * @brief Adjust perturbation cells so the net buoyant forcing is zero.
     * @param boxIdx Perturbation-box index on the level.
     * @param lev AMR level index.
     */
    void netZeroBuoyantAdjust (const int& boxIdx,
                               const int& lev)
    {
        // Creating local copy of PB box array and magnitude
        const amrex::BoxArray m_pb_ba = pb_ba[lev];
        for (amrex::MFIter mfi(pb_cell[lev], TileNoZ()) ; mfi.isValid(); ++mfi) {
            const amrex::Box& vbx = mfi.validbox();
            amrex::Box pbx = amrex::convert(m_pb_ba[boxIdx], vbx.ixType());
            amrex::Box ubx = pbx & vbx;
            if (ubx.ok()) {
               const amrex::Real adjust = tpi_pert_adjust;
               const amrex::Array4<amrex::Real>& pert_cell = pb_cell[lev].array(mfi);
               ParallelFor(ubx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
                   pert_cell(i,j,k) -= adjust;
               });
            }
        }
    }
 
// TODO: Test the difference between these two for Source term perturbation
#define USE_VOLUME_AVERAGE
    /**
     * @brief Compute mean horizontal velocity magnitude and direction for one perturbation box.
     * @param boxIdx Perturbation-box index on the level.
     * @param lev AMR level index.
     * @param mf_cons Conserved state used for iteration layout.
     * @param mf_xvel X-velocity field.
     * @param mf_yvel Y-velocity field.
     */
    void calc_tpi_meanMag_perBox (const int& boxIdx,
                                  const int& lev,
                                  amrex::MultiFab& mf_cons,
                                  amrex::MultiFab& mf_xvel,
                                  amrex::MultiFab& mf_yvel)
 
    {
        // Creating local copy of PB box array and magnitude
        const amrex::BoxArray m_pb_ba = pb_ba[lev];
        amrex::Real* m_pb_mag = pb_mag[lev].data();
        amrex::Real* m_pb_dir = pb_dir[lev].data();
        m_pb_mag[boxIdx] = zero; // Safety step
        m_pb_dir[boxIdx] = zero; // Safety step
 
        // Storage of averages per PB
        // Index: 0=u (vol/slab_lo), 1=v (vol/slab_lo)
        //        2=u (slab_hi),     3=v (slab_hi)
        int n_avg = 4;
        amrex::Vector<amrex::Real> avg_h(n_avg,zero);
        amrex::Gpu::DeviceVector<amrex::Real> avg_d(n_avg,zero);
        amrex::Real* avg = avg_d.data();
 
        // Averaging u & v components in single MFIter
        for (amrex::MFIter mfi(mf_cons, TileNoZ()); mfi.isValid(); ++mfi) {
 
            // CC valid box (inherited from mf_cons)
            const amrex::Box& vbx = mfi.validbox();
 
            // Box logic for u velocity
            auto ixtype_u    = mf_xvel.boxArray().ixType();
            amrex::Box vbx_u = amrex::convert(vbx,ixtype_u);
            amrex::Box pbx_u = amrex::convert(m_pb_ba[boxIdx], ixtype_u);
            amrex::Box ubx_u = pbx_u & vbx_u;
 
            // Box logic for v velocity
            auto ixtype_v    = mf_yvel.boxArray().ixType();
            amrex::Box vbx_v = amrex::convert(vbx,ixtype_v);
            amrex::Box pbx_v = amrex::convert(m_pb_ba[boxIdx], ixtype_v);
            amrex::Box ubx_v = pbx_v & vbx_v;
 
            // Operation over box union (U)
            if (ubx_u.ok()) {
                const amrex::Array4<const amrex::Real>& xvel_arry = mf_xvel.const_array(mfi);
 
                #ifdef USE_VOLUME_AVERAGE
                amrex::Real norm = one / static_cast<amrex::Real>(ubx_u.numPts());
                ParallelFor(amrex::Gpu::KernelInfo().setReduction(true), ubx_u, [=]
                AMREX_GPU_DEVICE(int i, int j, int k, amrex::Gpu::Handler const& handler) noexcept {
                    amrex::Gpu::deviceReduceSum(&avg[0], xvel_arry(i,j,k)*norm, handler);
                });
                #endif // USE_VOLUME_AVERAGE
 
                #ifdef USE_SLAB_AVERAGE
                amrex::Box ubxSlab_lo = makeSlab(ubx_u,2,ubx_u.smallEnd(2));
                amrex::Box ubxSlab_hi = makeSlab(ubx_u,2,ubx_u.bigEnd(2));
                amrex::Real norm_lo = one / static_cast<amrex::Real>(ubxSlab_lo.numPts());
                amrex::Real norm_hi = one / static_cast<amrex::Real>(ubxSlab_hi.numPts());
 
                // Average u in the low slab
                ParallelFor(amrex::Gpu::KernelInfo().setReduction(true), ubxSlab_lo, [=]
                AMREX_GPU_DEVICE(int i, int j, int k, amrex::Gpu::Handler const& handler) noexcept {
                    amrex::Gpu::deviceReduceSum(&avg[0], xvel_arry(i,j,k)*norm_lo, handler);
                });
 
                // Average u in the high slab
                ParallelFor(amrex::Gpu::KernelInfo().setReduction(true), ubxSlab_hi, [=]
                AMREX_GPU_DEVICE(int i, int j, int k, amrex::Gpu::Handler const& handler) noexcept {
                    amrex::Gpu::deviceReduceSum(&avg[2], xvel_arry(i,j,k)*norm_hi, handler);
                });
                #endif // USE_SLAB_AVERAGE
            } // if
 
            // Operation over box union (V)
            if (ubx_v.ok()) {
                const amrex::Array4<const amrex::Real>& yvel_arry = mf_yvel.const_array(mfi);
 
                #ifdef USE_VOLUME_AVERAGE
                amrex::Real norm = one / static_cast<amrex::Real>(ubx_v.numPts());
                ParallelFor(amrex::Gpu::KernelInfo().setReduction(true), ubx_v, [=]
                AMREX_GPU_DEVICE(int i, int j, int k, amrex::Gpu::Handler const& handler) noexcept {
                    amrex::Gpu::deviceReduceSum(&avg[1], yvel_arry(i,j,k)*norm, handler);
                });
                #endif // USE_VOLUME_AVERAGE
 
                #ifdef USE_SLAB_AVERAGE
                amrex::Box ubxSlab_lo = makeSlab(ubx_v,2,ubx_v.smallEnd(2));
                amrex::Box ubxSlab_hi = makeSlab(ubx_v,2,ubx_v.bigEnd(2));
                amrex::Real norm_lo = one / static_cast<amrex::Real>(ubxSlab_lo.numPts());
                amrex::Real norm_hi = one / static_cast<amrex::Real>(ubxSlab_hi.numPts());
 
                // Average v in the low slab
                ParallelFor(amrex::Gpu::KernelInfo().setReduction(true), ubxSlab_lo, [=]
                AMREX_GPU_DEVICE(int i, int j, int k, amrex::Gpu::Handler const& handler) noexcept {
                    amrex::Gpu::deviceReduceSum(&avg[1], yvel_arry(i,j,k)*norm_lo, handler);
                });
 
                // Average v in the high slab
                ParallelFor(amrex::Gpu::KernelInfo().setReduction(true), ubxSlab_hi, [=]
                AMREX_GPU_DEVICE(int i, int j, int k, amrex::Gpu::Handler const& handler) noexcept {
                    amrex::Gpu::deviceReduceSum(&avg[3], yvel_arry(i,j,k)*norm_hi, handler);
                });
                #endif // USE_SLAB_AVERAGE
            } // if
        } // MFIter
 
        // Copy from device back to host
        amrex::Gpu::copy(amrex::Gpu::deviceToHost, avg_d.begin(), avg_d.end(), avg_h.begin());
 
        // Computing the average magnitude within PB
        #ifdef USE_VOLUME_AVERAGE
        m_pb_mag[boxIdx] = std::sqrt(avg_h[0]*avg_h[0] + avg_h[1]*avg_h[1]);
        m_pb_dir[boxIdx] = std::atan(std::abs(avg_h[0]) / std::abs(avg_h[1]+std::numeric_limits<amrex::Real>::epsilon()));
        #endif
 
        #ifdef USE_SLAB_AVERAGE
        m_pb_mag[boxIdx] = myhalf*( std::sqrt(avg_h[0]*avg_h[0] + avg_h[1]*avg_h[1])
                                  + std::sqrt(avg_h[2]*avg_h[2] + avg_h[3]*avg_h[3]));
        m_pb_dir[boxIdx] = std::atan(std::abs(myhalf*(avg_h[0]+avg_h[2])) / std::abs(myhalf*(avg_h[1]+avg_h[3])+std::numeric_limits<amrex::Real>::epsilon()));
        #endif
    }
 
    /**
     * @brief Write perturbation debug information when debug output is enabled.
     */
    void debug (double /*time*/)
    {
        /*
        amrex::PrintToFile("BoxPerturbationOutput") << "#################### PB output at time = "
            << time << " ####################\n";
        amrex::PrintToFile("BoxPerturbationOutput") << " Using type: " <<  pt_type << "\n";
        amrex::PrintToFile("BoxPerturbationOutput") << " Net: " << tpi_net_buoyant << " Adjust : " <<  tpi_pert_adjust << "\n";
        for (int i = 0; i < pb_mag.size(); i++) {
            amrex::PrintToFile("BoxPerturbationOutput") << "[" << i
                << "] pb_Umag=" << pb_mag[i]
                << " | pb_interval=" << pb_interval[i]
                << " (" << pb_local_etime[i]
                << ") | pb_amp=" << pb_amp[i] << "\n";
        }
        amrex::PrintToFile("BoxPerturbationOutput") << "\n";
        */
    }
 
    amrex::Vector<int> pt_type;
 
    // Public data members
    amrex::Vector<amrex::BoxArray>                 pb_ba;  // PB box array
    amrex::Vector<amrex::Vector<amrex::Real>>      pb_mag; // BP mean magnitude [m/s]
    amrex::Vector<amrex::Vector<amrex::Real>>      pb_dir; // BP mean direction [deg]
 
    // Perturbation amplitude cell storage
    // This is after random assignment of equation (10) in Ma and Senocak 2023
    amrex::Vector<amrex::MultiFab>  pb_cell;
 
  private:
 
    // Private data members
    int tpi_layers;                   // Number of layers of perturbation boxes
    int tpi_offset;                   // Cells to offset the start of the perturbation region
 
    amrex::Vector<int> tpi_boxDim;    // Dimensions of each perturbation box
    amrex::Vector<int> tpi_direction; // Direction of the perturbation
 
    // Richardson Formulation
    amrex::Real tpi_nonDim;           // Richardson number
    amrex::Real tpi_Ti;               // Temperature intensity value
    amrex::Real tpi_Tinf;             // Reference temperature [K]
 
    // Physical dimensions
    amrex::Vector<amrex::Real> tpi_Hpb;              // PB height [m]
    amrex::Vector<amrex::Real> tpi_Lpb;              // PB length [m]
    amrex::Vector<amrex::Real> tpi_Wpb;              // PB width  [m]
    amrex::Vector<amrex::Real> tpi_lref;             // PB reference length [m]
 
    amrex::Real tpi_net_buoyant;      // Perturbation net-zero calculation storage
    amrex::Real tpi_pert_adjust;      // Perturbation adjust for net-zero per cell adjustment
 
    amrex::Vector<amrex::IntVect> ref_ratio;  // ref_ratio for multilevel run
    amrex::Real input_Ug;                     // input geostrophic wind speed to scale perturbations when using CPM
    amrex::Real input_w_amp;                  // input vertical velocity amplitude for w
    int         perturbation_klo;             // upper bound in Nz for perturbations
    int         perturbation_khi;             // lower bound in Nz for perturbations
 
    // Perturbation data arrays
    amrex::Vector<amrex::Vector<double>> pb_interval;         // PB update time [s]
    amrex::Vector<amrex::Vector<double>> pb_local_etime;      // PB local elapsed time [s]
    amrex::Vector<amrex::Vector<amrex::Real>> pb_amp;         // PB perturbation amplitude Ri:[K]
    amrex::Vector<amrex::Vector<amrex::Real>> pb_netZero;     // PB array used for net zero sum calculation
 
    /**
     * @brief Return a pseudo-random real value in a closed interval.
     * @param min Lower bound of the interval.
     * @param max Upper bound of the interval.
     * @return Pseudo-random value between min and max.
     */
    amrex::Real RandomReal (const amrex::Real min, const amrex::Real max)
    {
        amrex::Real r = (amrex::Real) rand() / (amrex::Real) RAND_MAX;
        return min + r * (max - min);
    }
 
};
 
/**
 * @var TurbulentPerturbation::pt_type
 * @brief Integer perturbation type for each AMR level.
 * @var TurbulentPerturbation::pb_ba
 * @brief Perturbation-box BoxArray for each AMR level.
 * @var TurbulentPerturbation::pb_mag
 * @brief Mean velocity magnitude for each perturbation box [m/s].
 * @var TurbulentPerturbation::pb_dir
 * @brief Mean velocity direction for each perturbation box.
 * @var TurbulentPerturbation::pb_cell
 * @brief Per-cell perturbation amplitude storage.
 * @var TurbulentPerturbation::tpi_layers
 * @brief Number of layers of perturbation boxes.
 * @var TurbulentPerturbation::tpi_offset
 * @brief Cell offset for the start of the perturbation region.
 * @var TurbulentPerturbation::tpi_boxDim
 * @brief Dimensions of each perturbation box.
 * @var TurbulentPerturbation::tpi_direction
 * @brief Boundary directions where perturbations are applied.
 * @var TurbulentPerturbation::tpi_nonDim
 * @brief Nondimensional number used by the perturbation formulation.
 * @var TurbulentPerturbation::tpi_Ti
 * @brief Temperature intensity used by the perturbation formulation.
 * @var TurbulentPerturbation::tpi_Tinf
 * @brief Reference temperature used by the perturbation formulation [K].
 * @var TurbulentPerturbation::tpi_Hpb
 * @brief Perturbation-box height for each AMR level [m].
 * @var TurbulentPerturbation::tpi_Lpb
 * @brief Perturbation-box length for each AMR level [m].
 * @var TurbulentPerturbation::tpi_Wpb
 * @brief Perturbation-box width for each AMR level [m].
 * @var TurbulentPerturbation::tpi_lref
 * @brief Perturbation-box reference length for each AMR level [m].
 * @var TurbulentPerturbation::tpi_net_buoyant
 * @brief Accumulated net buoyant perturbation used for correction.
 * @var TurbulentPerturbation::tpi_pert_adjust
 * @brief Per-cell perturbation adjustment used for net-zero buoyancy.
 * @var TurbulentPerturbation::ref_ratio
 * @brief Refinement ratios used by multilevel perturbation scaling.
 * @var TurbulentPerturbation::input_Ug
 * @brief Input geostrophic wind speed used to scale CPM perturbations.
 * @var TurbulentPerturbation::input_w_amp
 * @brief Input vertical velocity perturbation amplitude.
 * @var TurbulentPerturbation::perturbation_klo
 * @brief Lower vertical index bound for perturbations.
 * @var TurbulentPerturbation::perturbation_khi
 * @brief Upper vertical index bound for perturbations.
 * @var TurbulentPerturbation::pb_interval
 * @brief Perturbation update interval for each box [s].
 * @var TurbulentPerturbation::pb_local_etime
 * @brief Local elapsed time for each perturbation box [s].
 * @var TurbulentPerturbation::pb_amp
 * @brief Perturbation amplitude for each perturbation box.
 * @var TurbulentPerturbation::pb_netZero
 * @brief Per-box storage used for net-zero buoyancy calculation.
 */
#endif