TurbulentPerturbation Struct Reference

#include <ERF_TurbPertStruct.H>

Collaboration diagram for TurbulentPerturbation:

Public Member Functions
	~TurbulentPerturbation ()
void	init_tpi_type (const int lev, const PerturbationType &pert_type, const int max_level)
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)
void	calc_tpi_update (const int lev, const double dt, amrex::MultiFab &mf_xvel, amrex::MultiFab &mf_yvel, amrex::MultiFab &mf_cons)
void	apply_tpi (const int &lev, const amrex::Box &vbx, const int &comp, const amrex::IndexType &m_ixtype, const amrex::Array4< amrex::Real > &src_arr, const amrex::Array4< amrex::Real const > &pert_cell)
void	calc_tpi_amp (const int &lev, const int &boxIdx, const double &interval)
void	pseudoRandomPert (const int &boxIdx, const int &lev, const amrex::IndexType &m_ixtype)
void	zero_amp (const int &boxIdx, const int &lev, const amrex::IndexType &m_ixtype)
void	netZeroBuoyantAdd (const int &boxIdx, const int &lev)
void	netZeroBuoyantAdjust (const int &boxIdx, const int &lev)
void	calc_tpi_meanMag_perBox (const int &boxIdx, const int &lev, amrex::MultiFab &mf_cons, amrex::MultiFab &mf_xvel, amrex::MultiFab &mf_yvel)
void	debug (double)

Public Attributes
amrex::Vector< int >	pt_type
amrex::Vector< amrex::BoxArray >	pb_ba
amrex::Vector< amrex::Vector< amrex::Real > >	pb_mag
amrex::Vector< amrex::Vector< amrex::Real > >	pb_dir
amrex::Vector< amrex::MultiFab >	pb_cell

Private Member Functions
amrex::Real	RandomReal (const amrex::Real min, const amrex::Real max)

Private Attributes
int	tpi_layers
int	tpi_offset
amrex::Vector< int >	tpi_boxDim
amrex::Vector< int >	tpi_direction
amrex::Real	tpi_nonDim
amrex::Real	tpi_Ti
amrex::Real	tpi_Tinf
amrex::Vector< amrex::Real >	tpi_Hpb
amrex::Vector< amrex::Real >	tpi_Lpb
amrex::Vector< amrex::Real >	tpi_Wpb
amrex::Vector< amrex::Real >	tpi_lref
amrex::Real	tpi_net_buoyant
amrex::Real	tpi_pert_adjust
amrex::Vector< amrex::IntVect >	ref_ratio
amrex::Real	input_Ug
amrex::Real	input_w_amp
int	perturbation_klo
int	perturbation_khi
amrex::Vector< amrex::Vector< double > >	pb_interval
amrex::Vector< amrex::Vector< double > >	pb_local_etime
amrex::Vector< amrex::Vector< amrex::Real > >	pb_amp
amrex::Vector< amrex::Vector< amrex::Real > >	pb_netZero

Constructor & Destructor Documentation

~TurbulentPerturbation()

TurbulentPerturbation::~TurbulentPerturbation ( )

inline

{}

Member Function Documentation

apply_tpi()

void TurbulentPerturbation::apply_tpi ( const int &  lev, const amrex::Box &  vbx, const int &  comp, const amrex::IndexType &  m_ixtype, const amrex::Array4< amrex::Real > &  src_arr, const amrex::Array4< amrex::Real const > &  pert_cell  )

inline

    {
        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
                    });
                }
            }
        }
    }

Referenced by make_sources().

Here is the call graph for this function:

Here is the caller graph for this function:

calc_tpi_amp()

void TurbulentPerturbation::calc_tpi_amp ( const int &  lev, const int &  boxIdx, const double &  interval  )

inline

    {
        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);
            }
        }
    }

Referenced by calc_tpi_update().

Here is the caller graph for this function:

calc_tpi_meanMag_perBox()

void TurbulentPerturbation::calc_tpi_meanMag_perBox ( const int &  boxIdx, const int &  lev, amrex::MultiFab &  mf_cons, amrex::MultiFab &  mf_xvel, amrex::MultiFab &  mf_yvel  )

inline

    {
        // 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
    }

Referenced by calc_tpi_update().

Here is the call graph for this function:

Here is the caller graph for this function:

calc_tpi_update()

void TurbulentPerturbation::calc_tpi_update ( const int  lev, const double  dt, amrex::MultiFab &  mf_xvel, amrex::MultiFab &  mf_yvel, amrex::MultiFab &  mf_cons  )

inline

    {
        // 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);
            }
        }
    }

Here is the call graph for this function:

debug()

void TurbulentPerturbation::debug ( double  )

inline

    {
        /*
        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";
        */
    }

init_tpi()

void TurbulentPerturbation::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  )

inline

    {
        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";
        */
    }

Here is the call graph for this function:

init_tpi_type()

void TurbulentPerturbation::init_tpi_type ( const int  lev, const PerturbationType &  pert_type, const int  max_level  )

inline

    {
        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;
        }
    }

netZeroBuoyantAdd()

void TurbulentPerturbation::netZeroBuoyantAdd ( const int &  boxIdx, const int &  lev  )

inline

    {
        // 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];
            }
        }
    }

Referenced by calc_tpi_update().

Here is the call graph for this function:

Here is the caller graph for this function:

netZeroBuoyantAdjust()

void TurbulentPerturbation::netZeroBuoyantAdjust ( const int &  boxIdx, const int &  lev  )

inline

    {
        // 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;
               });
            }
        }
    }

Referenced by calc_tpi_update().

Here is the call graph for this function:

Here is the caller graph for this function:

pseudoRandomPert()

void TurbulentPerturbation::pseudoRandomPert ( const int &  boxIdx, const int &  lev, const amrex::IndexType &  m_ixtype  )

inline

    {
        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;
                    });
                }
            }
        }
    }

Referenced by calc_tpi_update().

Here is the call graph for this function:

Here is the caller graph for this function:

RandomReal()

amrex::Real TurbulentPerturbation::RandomReal ( const amrex::Real  min, const amrex::Real  max  )

inlineprivate

    {
        amrex::Real r = (amrex::Real) rand() / (amrex::Real) RAND_MAX;
        return min + r * (max - min);
    }

Referenced by calc_tpi_update(), and pseudoRandomPert().

Here is the caller graph for this function:

zero_amp()

void TurbulentPerturbation::zero_amp ( const int &  boxIdx, const int &  lev, const amrex::IndexType &  m_ixtype  )

inline

    {
 
        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;
                });
            }
        }
    }

Referenced by calc_tpi_update().

Here is the call graph for this function:

Here is the caller graph for this function:

Member Data Documentation

input_Ug

amrex::Real TurbulentPerturbation::input_Ug

private

Referenced by calc_tpi_amp(), and init_tpi().

input_w_amp

amrex::Real TurbulentPerturbation::input_w_amp

private

Referenced by calc_tpi_amp(), and init_tpi().

pb_amp

amrex::Vector<amrex::Vector<amrex::Real> > TurbulentPerturbation::pb_amp

private

Referenced by calc_tpi_amp(), init_tpi(), and pseudoRandomPert().

pb_ba

amrex::Vector<amrex::BoxArray> TurbulentPerturbation::pb_ba

Referenced by apply_tpi(), calc_tpi_meanMag_perBox(), calc_tpi_update(), init_tpi(), netZeroBuoyantAdd(), netZeroBuoyantAdjust(), pseudoRandomPert(), and zero_amp().

pb_cell

amrex::Vector<amrex::MultiFab> TurbulentPerturbation::pb_cell

Referenced by init_tpi(), make_sources(), netZeroBuoyantAdd(), netZeroBuoyantAdjust(), pseudoRandomPert(), and zero_amp().

pb_dir

amrex::Vector<amrex::Vector<amrex::Real> > TurbulentPerturbation::pb_dir

Referenced by calc_tpi_meanMag_perBox(), calc_tpi_update(), and init_tpi().

pb_interval

amrex::Vector<amrex::Vector<double> > TurbulentPerturbation::pb_interval

private

Referenced by calc_tpi_update(), and init_tpi().

pb_local_etime

amrex::Vector<amrex::Vector<double> > TurbulentPerturbation::pb_local_etime

private

Referenced by calc_tpi_update(), and init_tpi().

pb_mag

amrex::Vector<amrex::Vector<amrex::Real> > TurbulentPerturbation::pb_mag

Referenced by calc_tpi_amp(), calc_tpi_meanMag_perBox(), calc_tpi_update(), and init_tpi().

pb_netZero

amrex::Vector<amrex::Vector<amrex::Real> > TurbulentPerturbation::pb_netZero

private

Referenced by calc_tpi_update(), init_tpi(), and netZeroBuoyantAdd().

perturbation_khi

int TurbulentPerturbation::perturbation_khi

private

Referenced by init_tpi().

perturbation_klo

int TurbulentPerturbation::perturbation_klo

private

Referenced by init_tpi().

pt_type

amrex::Vector<int> TurbulentPerturbation::pt_type

Referenced by apply_tpi(), calc_tpi_amp(), calc_tpi_update(), init_tpi(), init_tpi_type(), pseudoRandomPert(), and zero_amp().

ref_ratio

amrex::Vector<amrex::IntVect> TurbulentPerturbation::ref_ratio

private

Referenced by calc_tpi_amp(), and init_tpi().

tpi_boxDim

amrex::Vector<int> TurbulentPerturbation::tpi_boxDim

private

Referenced by init_tpi().

tpi_direction

amrex::Vector<int> TurbulentPerturbation::tpi_direction

private

Referenced by init_tpi().

tpi_Hpb

amrex::Vector<amrex::Real> TurbulentPerturbation::tpi_Hpb

private

Referenced by calc_tpi_amp(), and init_tpi().

tpi_layers

int TurbulentPerturbation::tpi_layers

private

Referenced by calc_tpi_update(), and init_tpi().

tpi_Lpb

amrex::Vector<amrex::Real> TurbulentPerturbation::tpi_Lpb

private

Referenced by calc_tpi_update(), and init_tpi().

tpi_lref

amrex::Vector<amrex::Real> TurbulentPerturbation::tpi_lref

private

Referenced by calc_tpi_update(), and init_tpi().

tpi_net_buoyant

amrex::Real TurbulentPerturbation::tpi_net_buoyant

private

Referenced by calc_tpi_update(), and init_tpi().

tpi_nonDim

amrex::Real TurbulentPerturbation::tpi_nonDim

private

Referenced by calc_tpi_amp(), and init_tpi().

tpi_offset

int TurbulentPerturbation::tpi_offset

private

Referenced by init_tpi().

tpi_pert_adjust

amrex::Real TurbulentPerturbation::tpi_pert_adjust

private

Referenced by calc_tpi_update(), init_tpi(), and netZeroBuoyantAdjust().

tpi_Ti

amrex::Real TurbulentPerturbation::tpi_Ti

private

Referenced by calc_tpi_amp(), and init_tpi().

tpi_Tinf

amrex::Real TurbulentPerturbation::tpi_Tinf

private

Referenced by calc_tpi_amp(), and init_tpi().

tpi_Wpb

amrex::Vector<amrex::Real> TurbulentPerturbation::tpi_Wpb

private

Referenced by calc_tpi_update(), and init_tpi().

---

The documentation for this struct was generated from the following file:

• Source/DataStructs/ERF_TurbPertStruct.H