Runtime state and operations for turbulent perturbation forcing. More...

#include <ERF_TurbPertStruct.H>

Collaboration diagram for TurbulentPerturbation:

[legend]

Classes
struct	UniformClock

Public Member Functions
	~TurbulentPerturbation ()
	Destroy turbulent perturbation state. More...

void	init_tpi_type (const int lev, const PerturbationType &pert_type, const int max_level)
	Store the integer perturbation-type selector for one AMR level. More...

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, const bool automatic=false)
	Initialize turbulent perturbation regions and per-box storage. More...

amrex::Box	automatic_face_box (const amrex::Box &subdomain, int face) const

bool	automatic_update_due (int lev, double time) const

bool	uniform_refresh_enabled (int lev) const

bool	uniform_wind_update_due (int lev, double time) const

amrex::Real	uniform_sampling_height () const

amrex::Real	uniform_sampling_period () const

bool	uniform_sampling_face (int face) const

void	set_uniform_advection (int lev, double time, amrex::Real speed, double wind_period)

double	uniform_refresh_interval (int lev) const

bool	calc_tpi_uniform_update (int lev, double time)

void	configure_automatic (int lev, double time, double period, const amrex::Vector< AutomaticPerturbationFace > &faces, const amrex::Geometry &geom, const amrex::MultiFab z_cc, const amrex::MultiFab z_nd, const amrex::MultiFab *z_ground)

void	calc_tpi_update (const int lev, const double dt, amrex::MultiFab &mf_xvel, amrex::MultiFab &mf_yvel, amrex::MultiFab &mf_cons)
	Update perturbation amplitudes and intervals when each box is ready. More...

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, const amrex::Real scale=one)
	Apply stored turbulent perturbations to a source or state array. More...

void	calc_tpi_amp (const int &lev, const int &boxIdx, const double &interval)
	Compute the perturbation amplitude for one perturbation box. More...

void	pseudoRandomPert (const int &boxIdx, const int &lev, const amrex::IndexType &m_ixtype)
	Assign pseudo-random perturbations to cells in one perturbation box. More...

void	zero_amp (const int &boxIdx, const int &lev, const amrex::IndexType &m_ixtype)
	Reset CPM perturbation amplitudes in one perturbation box to zero. More...

void	netZeroBuoyant (const int &lev)
	Remove the mean of the box-perturbation field, so the net buoyant forcing is zero. More...

void	calc_tpi_meanMag (const int &lev, const amrex::Vector< int > &update_box, amrex::MultiFab &mf_cons, amrex::MultiFab &mf_xvel, amrex::MultiFab &mf_yvel)
	Compute the mean horizontal velocity magnitude and direction of the refreshing perturbation boxes. More...

void	debug (double)
	Write perturbation debug information when debug output is enabled. More...

Static Public Member Functions
static amrex::Real	box_sum (const amrex::Array4< const amrex::Real > &arr, const amrex::Box &bx)
	Sum an array over a box. More...

Public Attributes
amrex::Vector< int >	pt_type
	Integer perturbation type for each AMR level. More...

amrex::Vector< amrex::BoxArray >	pb_ba
	Perturbation-box BoxArray for each AMR level. More...

amrex::Vector< amrex::Vector< amrex::Real > >	pb_mag
	Mean velocity magnitude for each perturbation box [m/s]. More...

amrex::Vector< amrex::Vector< amrex::Real > >	pb_dir
	Mean velocity direction for each perturbation box. More...

amrex::Vector< amrex::MultiFab >	pb_cell
	Per-cell perturbation amplitude storage. More...

Private Attributes
bool	uniform_mode {true}

amrex::Real	uniform_height {amrex::Real(20)}

amrex::Real	uniform_wind_period {amrex::Real(600)}

amrex::Vector< UniformClock >	uniform_clock

bool	auto_mode {false}

amrex::Vector< amrex::Vector< amrex::Real > >	pb_Ug

amrex::Vector< amrex::MultiFab >	auto_mask

amrex::Vector< double >	auto_next_update

amrex::Vector< int >	auto_ready

int	tpi_layers
	Number of layers of perturbation boxes. More...

int	tpi_offset
	Cell offset for the start of the perturbation region. More...

amrex::Vector< int >	tpi_boxDim
	Dimensions of each perturbation box. More...

amrex::Vector< int >	tpi_direction
	Boundary directions where perturbations are applied. More...

amrex::Real	tpi_nonDim
	Nondimensional number used by the perturbation formulation. More...

amrex::Real	tpi_Ti
	Temperature intensity used by the perturbation formulation. More...

amrex::Real	tpi_Tinf
	Reference temperature used by the perturbation formulation [K]. More...

amrex::Vector< amrex::Real >	tpi_Hpb
	Perturbation-box height for each AMR level [m]. More...

amrex::Vector< amrex::Real >	tpi_Lpb
	Perturbation-box length for each AMR level [m]. More...

amrex::Vector< amrex::Real >	tpi_Wpb
	Perturbation-box width for each AMR level [m]. More...

amrex::Vector< amrex::Real >	tpi_lref
	Perturbation-box reference length for each AMR level [m]. More...

amrex::Real	tpi_net_buoyant
	Accumulated net buoyant perturbation used for correction. More...

amrex::Real	tpi_pert_adjust
	Per-cell perturbation adjustment used for net-zero buoyancy. More...

amrex::Vector< amrex::IntVect >	ref_ratio
	Refinement ratios used by multilevel perturbation scaling. More...

amrex::Real	input_Ug
	Input geostrophic wind speed used to scale CPM perturbations. More...

amrex::Real	input_w_amp
	Input vertical velocity perturbation amplitude. More...

int	perturbation_klo
	Lower vertical index bound for perturbations. More...

int	perturbation_khi
	Upper vertical index bound for perturbations. More...

amrex::Vector< amrex::Vector< double > >	pb_interval
	Perturbation update interval for each box [s]. More...

amrex::Vector< amrex::Vector< double > >	pb_local_etime
	Local elapsed time for each perturbation box [s]. More...

amrex::Vector< amrex::Vector< amrex::Real > >	pb_amp
	Perturbation amplitude for each perturbation box. More...

amrex::Vector< amrex::Vector< amrex::Real > >	pb_netZero
	Per-box storage used for net-zero buoyancy calculation. More...

amrex::Vector< amrex::Vector< int > >	pb_nupdate
	Number of times each perturbation box has been refreshed; keys its random numbers. More...

std::uint64_t	tpi_seed = 0
	Seed of the hashed random numbers (1024 with erf.fix_random_seed, else erf.random_seed, else 0). More...

Static Private Attributes
static constexpr int	hash_comp_cell = 0
	Hash salt of the per-cell perturbation values. More...

static constexpr int	hash_comp_cpm = 1
	Hash salt of the per-box CPM values. More...

static constexpr int	hash_comp_interval = 2
	Hash salt of the refresh-interval variation. More...

Detailed Description

Runtime state and operations for turbulent perturbation forcing.

Constructor & Destructor Documentation

◆ ~TurbulentPerturbation()

TurbulentPerturbation::~TurbulentPerturbation ( )

inline

Destroy turbulent perturbation state.

49 {}

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,
		const amrex::Real	scale = `one`
	)

inline

Apply stored turbulent perturbations to a source or state array.

Parameters

lev	AMR level index.
vbx	Valid box over which perturbations may be applied.
comp	Component to modify, or -1 for vertical velocity perturbations.
m_ixtype	Index type of the destination array.
src_arr	Destination array receiving perturbation increments.
pert_cell	Stored perturbation amplitudes.
scale	Factor applied to the stored amplitudes (the time step for the direct perturbation, whose stored amplitudes are rates).

     {
         if (uniform_refresh_enabled(lev) && !uniform_clock[lev].apply) { return; }
         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) += scale * pert_cell(i,j,k);
                     });
                 } else {
                     ParallelFor(ubx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
                         src_arr(i,j,k,comp) += scale * 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:

◆ automatic_face_box()

amrex::Box TurbulentPerturbation::automatic_face_box	(	const amrex::Box &	subdomain,
		int	face
	)		const

inline

     {
         amrex::Box region = subdomain;
         region.grow(0, -tpi_offset);
         region.grow(1, -tpi_offset);
         const int dir = face % 2;
         const int width = tpi_layers * tpi_boxDim[dir];
         if (face < 2) {
             region.setBig(dir, std::min(region.bigEnd(dir), region.smallEnd(dir)+width-1));
         } else {
             region.setSmall(dir, std::max(region.smallEnd(dir), region.bigEnd(dir)-width+1));
         }
         return region;
     }

Referenced by init_tpi().

Here is the caller graph for this function:

◆ automatic_update_due()

bool TurbulentPerturbation::automatic_update_due	(	int	lev,
		double	time
	)		const

inline

     {
         return !auto_ready[lev] || time >= auto_next_update[lev];
     }

◆ box_sum()

static amrex::Real TurbulentPerturbation::box_sum	(	const amrex::Array4< const amrex::Real > &	arr,
		const amrex::Box &	bx
	)

inlinestatic

Sum an array over a box.

Parameters

arr	Array to sum.
bx	Box to sum over.

Returns: Sum of arr over bx.

     {
         amrex::ReduceOps<amrex::ReduceOpSum> reduce_op;
         amrex::ReduceData<amrex::Real> reduce_data(reduce_op);
         using ReduceTuple = typename decltype(reduce_data)::Type;
         reduce_op.eval(bx, reduce_data, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept -> ReduceTuple
         {
             return {arr(i,j,k)};
         });
         return amrex::get<0>(reduce_data.value(reduce_op));
     }

Referenced by calc_tpi_meanMag(), and netZeroBuoyant().

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

Compute the perturbation amplitude for one perturbation box.

Parameters

lev	AMR level index.
boxIdx	Perturbation-box index on the level.
interval	Current perturbation update interval [s].

     {
         pb_amp[lev][boxIdx] = zero; // Safety step
         if (pt_type[lev] == 2) { // CPM
             constexpr amrex::Real cp = amrex::Real(1004); // specific heat of air [J/(kg K)]
             constexpr amrex::Real Ec = amrex::Real(0.2); // Eckert number
             const amrex::Real Ug = auto_mode ? pb_Ug[lev][boxIdx] : input_Ug;
             pb_amp[lev][boxIdx] = (Ug * 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];
  
             // No mean wind: nothing to trip, and with tpi_Ti > 0 the amplitude below would be
             //    0/0. The box is not refreshed again, so a zero amplitude also clears the
             //    perturbation it held
             if (Um == zero) { return; }
  
             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;
  
             // Divide by the interval to get a forcing rate [K/s], Ma & Senocak (2023) Eq. 7.
             //    The source perturbation adds this rate as a source term; the direct
             //    perturbation adds it times the time step to the state at every step. Either
             //    way the total over one refresh interval is delta phi, even if dt changes
             pb_amp[lev][boxIdx] /= static_cast<amrex::Real>(interval);
         }
     }

Referenced by calc_tpi_uniform_update(), and calc_tpi_update().

Here is the caller graph for this function:

◆ calc_tpi_meanMag()

void TurbulentPerturbation::calc_tpi_meanMag	(	const int &	lev,
		const amrex::Vector< int > &	update_box,
		amrex::MultiFab &	mf_cons,
		amrex::MultiFab &	mf_xvel,
		amrex::MultiFab &	mf_yvel
	)

inline

Compute the mean horizontal velocity magnitude and direction of the refreshing perturbation boxes.

The velocities are summed over every grid box a perturbation box overlaps, each face once, and over all ranks, so every rank gets the same mean whatever the domain decomposition.

Parameters

lev	AMR level index.
update_box	Nonzero for the perturbation boxes that refresh this step.
mf_cons	Conserved state used for iteration layout.
mf_xvel	X-velocity field.
mf_yvel	Y-velocity field.

     {
         const amrex::BoxArray& m_pb_ba = pb_ba[lev];
         const int nboxes = m_pb_ba.size();
  
         // Per box: sum of u, number of u faces, sum of v, number of v faces
         amrex::Vector<amrex::Real> sums(4*nboxes, zero);
         for (amrex::MFIter mfi(mf_cons); mfi.isValid(); ++mfi) {
             const amrex::Box& vbx = mfi.validbox(); // cell centered
             const amrex::Array4<const amrex::Real>& xvel_arry = mf_xvel.const_array(mfi);
             const amrex::Array4<const amrex::Real>& yvel_arry = mf_yvel.const_array(mfi);
             for (const auto& isect : m_pb_ba.intersections(vbx)) {
                 const int boxIdx = isect.first;
                 if (!update_box[boxIdx]) { continue; }
                 const amrex::Box& pbx = m_pb_ba[boxIdx];
  
                 // The faces of the overlap's cells; the high face only where the overlap ends
                 // at the perturbation box's high side, so a face shared by two grid boxes counts once
                 amrex::Box ubx_u = amrex::surroundingNodes(isect.second, 0);
                 if (isect.second.bigEnd(0) < pbx.bigEnd(0)) { ubx_u.growHi(0, -1); }
                 amrex::Box ubx_v = amrex::surroundingNodes(isect.second, 1);
                 if (isect.second.bigEnd(1) < pbx.bigEnd(1)) { ubx_v.growHi(1, -1); }
  
                 sums[4*boxIdx  ] += box_sum(xvel_arry, ubx_u);
                 sums[4*boxIdx+1] += static_cast<amrex::Real>(ubx_u.numPts());
                 sums[4*boxIdx+2] += box_sum(yvel_arry, ubx_v);
                 sums[4*boxIdx+3] += static_cast<amrex::Real>(ubx_v.numPts());
             }
         }
         amrex::ParallelAllReduce::Sum(sums.data(), static_cast<int>(sums.size()),
                                       amrex::ParallelContext::CommunicatorSub());
  
         for (int boxIdx = 0; boxIdx < nboxes; boxIdx++) {
             if (!update_box[boxIdx]) { continue; }
             const amrex::Real u_avg = (sums[4*boxIdx+1] > zero) ? sums[4*boxIdx  ] / sums[4*boxIdx+1] : zero;
             const amrex::Real v_avg = (sums[4*boxIdx+3] > zero) ? sums[4*boxIdx+2] / sums[4*boxIdx+3] : zero;
             pb_mag[lev][boxIdx] = std::sqrt(u_avg*u_avg + v_avg*v_avg);
             pb_dir[lev][boxIdx] = std::atan(std::abs(u_avg) / std::abs(v_avg+std::numeric_limits<amrex::Real>::epsilon()));
         }
     }

Referenced by calc_tpi_update().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ calc_tpi_uniform_update()

bool TurbulentPerturbation::calc_tpi_uniform_update	(	int	lev,
		double	time
	)

inline

     {
         AMREX_ALWAYS_ASSERT(uniform_refresh_enabled(lev));
         auto& clock = uniform_clock[lev];
         clock.apply = false;
         if ((auto_mode && !auto_ready[lev]) || clock.interval <= 0.0 ||
             time < clock.next_refresh || pb_ba[lev].empty()) { return false; }
  
         clock.next_refresh = time + clock.interval;
         clock.apply = true;
         pb_cell[lev].setVal(zero);
         const auto m_ixtype = pb_cell[lev].ixType();
         for (int b = 0; b < pb_ba[lev].size(); ++b) {
             pb_interval[lev][b] = clock.interval / double(tpi_layers);
             pb_local_etime[lev][b] = zero;
             calc_tpi_amp(lev, b, pb_interval[lev][b]);
             pseudoRandomPert(b, lev, m_ixtype);
             ++pb_nupdate[lev][b];
         }
         return true;
     }

Here is the call 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

Update perturbation amplitudes and intervals when each box is ready.

Parameters

lev	AMR level index.
dt	Time step for elapsed-time accumulation.
mf_xvel	X-velocity field used to compute box mean velocities.
mf_yvel	Y-velocity field used to compute box mean velocities.
mf_cons	Conserved state used for iteration layout.

     {
         auto m_ixtype = mf_cons.boxArray().ixType(); // safety step
         if (pt_type[lev] == 3) { m_ixtype = amrex::IndexType(amrex::IntVect(0,0,1)); }
  
         const int nboxes = pb_ba[lev].size();
         if (nboxes == 0 || (auto_mode && !auto_ready[lev])) { return; }
  
         // Which boxes refresh this step. pb_interval and pb_local_etime hold the same values on
         // every rank, so every rank makes the same choice.
         amrex::Vector<int> update_box(nboxes, 0);
         for (int boxIdx = 0; boxIdx < nboxes; boxIdx++) {
             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[boxIdx] = ((pb_local_etime[lev][boxIdx] >= pb_interval[lev][boxIdx]*tpi_layers) && pb_local_etime[lev][boxIdx] != -one);
             } else {
                 update_box[boxIdx] = ( pb_local_etime[lev][boxIdx] >= pb_interval[lev][boxIdx] && pb_local_etime[lev][boxIdx] != -one );
             }
         }
  
         // Compute mean velocity of each refreshing perturbation box
         if (!auto_mode) { calc_tpi_meanMag(lev, update_box, mf_cons, mf_xvel, mf_yvel); }
  
         for (int boxIdx = 0; boxIdx < nboxes; boxIdx++) {
             if ( update_box[boxIdx] ) {
  
                 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];
                     }
                     // 10% variation
                     const amrex::Real r = erf_hash_rng::hash_uniform(boxIdx, pb_nupdate[lev][boxIdx], 0,
                                                                      hash_comp_interval, lev, tpi_seed);
                     pb_interval[lev][boxIdx] = (amrex::Real(0.9) + amrex::Real(0.2)*r) * interval;
  
                     // Reset local elapsed time
                     pb_local_etime[lev][boxIdx] = zero;
                 } else {
                     // No mean wind through this box, so no advective time scale: it is not refreshed again
                     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);
  
                 pb_nupdate[lev][boxIdx]++;
  
             } 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
                 if (pb_local_etime[lev][boxIdx] != -one) { pb_local_etime[lev][boxIdx] += dt; }
             } // if
         } // for
  
         if (pt_type[lev] < 2) { // box perturbation method only
             netZeroBuoyant(lev);
         }
     }

Here is the call graph for this function:

◆ configure_automatic()

void TurbulentPerturbation::configure_automatic	(	int	lev,
		double	time,
		double	period,
		const amrex::Vector< AutomaticPerturbationFace > &	faces,
		const amrex::Geometry &	geom,
		const amrex::MultiFab *	z_cc,
		const amrex::MultiFab *	z_nd,
		const amrex::MultiFab *	z_ground
	)

inline

     {
         auto_mask[lev].setVal(zero);
         pb_cell[lev].setVal(zero);
         pb_mag[lev].assign(pb_ba[lev].size(), zero);
         pb_Ug[lev].assign(pb_ba[lev].size(), zero);
         pb_dir[lev].assign(pb_ba[lev].size(), zero);
         const bool on_w_faces = (pt_type[lev] == 3);
         const bool terrain = (z_cc != nullptr);
         AMREX_ALWAYS_ASSERT(!terrain || (z_nd && z_ground));
         const int klo = geom.Domain().smallEnd(2);
         const auto dx = geom.CellSizeArray();
         const auto problo = geom.ProbLoArray();
         for (const auto& face : faces) {
             if (!face.inflow()) { continue; }
             const amrex::Real speed = std::hypot(face.u, face.v);
             for (const auto& isect : pb_ba[lev].intersections(face.region)) {
                 const int b = isect.first;
                 pb_Ug[lev][b] = std::max(pb_Ug[lev][b], face.Ug);
                 if (speed > pb_mag[lev][b]) {
                     pb_mag[lev][b] = speed;
                     pb_dir[lev][b] = std::atan2(std::abs(face.u), std::abs(face.v));
                 }
             }
             const amrex::Real height = face.height;
             const amrex::Box region = amrex::convert(face.region, auto_mask[lev].ixType());
             for (amrex::MFIter mfi(auto_mask[lev]); mfi.isValid(); ++mfi) {
                 const amrex::Box bx = mfi.validbox() & region;
                 if (!bx.ok()) { continue; }
                 const auto mask = auto_mask[lev].array(mfi);
                 const auto cc = terrain ? z_cc->const_array(mfi) : amrex::Array4<const amrex::Real>{};
                 const auto nd = terrain ? z_nd->const_array(mfi) : amrex::Array4<const amrex::Real>{};
                 const auto ground = terrain ? z_ground->const_array(mfi) : amrex::Array4<const amrex::Real>{};
                 ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
                     amrex::Real z = problo[2] + (k-klo + (on_w_faces ? zero : myhalf))*dx[2];
                     amrex::Real zs = problo[2];
                     if (terrain) {
                         zs = fourth*(ground(i,j,klo) + ground(i+1,j,klo) +
                                      ground(i,j+1,klo) + ground(i+1,j+1,klo));
                         z = on_w_faces ? fourth*(nd(i,j,k) + nd(i+1,j,k) +
                                                  nd(i,j+1,k) + nd(i+1,j+1,k))
                                        : cc(i,j,k);
                     }
                     if (z > zs && z-zs <= height) { mask(i,j,k) = one; }
                 });
             }
         }
         // Reconsider even boxes that previously had calm winds. They must be able
         // to start again after a later parent update changes the inflow.
         pb_interval[lev].assign(pb_ba[lev].size(), -one);
         pb_local_etime[lev].assign(pb_ba[lev].size(), zero);
         auto_ready[lev] = 1;
         auto_next_update[lev] = time + period;
     }

Here is the call graph for this function:

◆ debug()

void TurbulentPerturbation::debug ( double )

inline

Write perturbation debug information when debug output is enabled.

     {
         /*
         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,
		const bool	automatic = `false`
	)

inline

Initialize turbulent perturbation regions and per-box storage.

Parameters

lev	AMR level index.
subdomains_lev	Rectangular level subdomains used to define perturbation regions.
dx	Cell spacing for this level.
ba	State BoxArray used for perturbation-cell storage.
dm	DistributionMapping for perturbation-cell storage.
ngrow_state	Number of ghost cells for perturbation-cell storage.
pp_prefix	ParmParse prefix for the ERF input namespace.
refRatio	Refinement ratios between AMR levels.
max_level	Maximum AMR level configured for the run.
automatic	Select runtime inflow/height masks driven by the parent.

     {
         ref_ratio = refRatio;
         auto_mode = automatic;
  
         amrex::ParmParse pp(pp_prefix);
  
         uniform_mode = true;
         pp.queryAdd("perturbation_uniform_refresh", uniform_mode);
         uniform_height = amrex::Real(20);
         pp.queryAdd("perturbation_refresh_height", uniform_height);
         uniform_wind_period = amrex::Real(600);
         pp.queryAdd("perturbation_refresh_wind_period", uniform_wind_period);
         if (uniform_refresh_enabled(lev) &&
             (!(uniform_height > zero) || !std::isfinite(uniform_height) ||
              !(uniform_wind_period > zero) || !std::isfinite(uniform_wind_period))) {
             amrex::Abort("Uniform CPM refresh requires positive, finite "
                          "perturbation_refresh_height and perturbation_refresh_wind_period");
         }
  
         // Reading inputs, and placing assertion for the perturbation inflow to work
         pp.getarr("perturbation_box_dims",tpi_boxDim);
         if (auto_mode) {
             tpi_direction = {1, 1, 0, 1, 1, 0};
         } else {
             pp.getarr("perturbation_direction",tpi_direction);
         }
         if (tpi_boxDim.size() != 3) {
             amrex::Abort("perturbation_box_dims must have 3 entries (Nx Ny Nz)");
         }
         if (tpi_direction.size() != 6) {
             amrex::Abort("perturbation_direction must have 6 entries (west south bottom east north top)");
         }
         pp.get("perturbation_layers",tpi_layers);
         pp.get("perturbation_offset",tpi_offset);
  
         tpi_nonDim = zero;
         pp.queryAdd("perturbation_nondimensional",tpi_nonDim);
  
         tpi_Tinf = amrex::Real(300.);
         pp.queryAdd("perturbation_T_infinity",tpi_Tinf);
  
         tpi_Ti = zero;
         pp.queryAdd("perturbation_T_intensity",tpi_Ti);
  
         input_Ug = zero;
         pp.queryAdd("perturbation_Ug",input_Ug);
  
         input_w_amp = zero;
         pp.queryAdd("perturbation_w_amp",input_w_amp);
  
         // Seed of the perturbation's random numbers, read the way erf.deterministic_ic_pert reads
         // it. The numbers are hashed from the global cell index (or the perturbation-box index) and
         // the box's refresh count, so they do not depend on the domain decomposition (#4317).
         {
             int fix_random_seed = 0;
             long random_seed = -1;
             amrex::ParmParse pp_erf("erf");
             pp_erf.query("fix_random_seed", fix_random_seed);
             pp_erf.query("random_seed", random_seed);
             tpi_seed = fix_random_seed ? std::uint64_t(1024)
                      : ((random_seed >= 0) ? static_cast<std::uint64_t>(random_seed) : std::uint64_t(0));
         }
  
         // Negative sentinels rather than the queryAdd return values: k = 0 is a legal
         // bound, so the natural default cannot double as "unset", and the return value
         // only reports whether the key existed before the call.  init_tpi runs once per
         // level, so testing the return would make every level above 0 see "specified"
         // with the value level 0 planted, skipping the valid_box_lo/hi defaults below.
         perturbation_klo = -1;
         pp.queryAdd("perturbation_klo", perturbation_klo);
         const bool have_klo = (perturbation_klo >= 0);
  
         perturbation_khi = -1;
         pp.queryAdd("perturbation_khi", perturbation_khi);
         const bool have_khi = (perturbation_khi >= 0);
  
         // 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 ((auto_mode || uniform_refresh_enabled(lev)) && tpi_layers == 0) {
             amrex::Abort("Automatic perturbations and uniform CPM refresh require perturbation_layers > 0");
         }
         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();
  
             if (auto_mode) {
                 // Allocate all possible edges once. A mask chooses inflow and height
                 // at runtime, including corners when either adjoining edge is inflow.
                 for (int face = 0; face < 4; ++face) {
                     const amrex::Box region = automatic_face_box(subdomain_box, face);
                     if (region.ok()) { tmp_bl.push_back(region); }
                 }
                 continue;
             }
  
             // 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);
         if (auto_mode) { tmp_ba.removeOverlap(); }
         tmp_ba.maxSize(boxSize);
         if (auto_mode && tmp_ba.empty()) {
             amrex::Abort("Automatic perturbation strips are empty; reduce perturbation_offset");
         }
  
         const int num_levels = max_level + 1;
         if (pb_ba.size() < num_levels) {
             pb_ba.resize(num_levels);
             pb_mag.resize(num_levels);
             pb_Ug.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_nupdate.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);
             auto_mask.resize(num_levels);
             auto_next_update.resize(num_levels, -1.0);
             auto_ready.resize(num_levels, 0);
             uniform_clock.resize(num_levels);
         }
         uniform_clock[lev] = UniformClock{};
  
         pb_ba[lev] = tmp_ba;
  
         // Initializing mean magnitude and direction vectors
         pb_mag[lev].resize(pb_ba[lev].size(), zero);
         pb_Ug[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);
         pb_nupdate[lev].resize(pb_ba[lev].size(), 0);
  
         // 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.);
         if (auto_mode) {
             auto_mask[lev].define(pb_cell[lev].boxArray(), dm, 1, 0);
             auto_mask[lev].setVal(zero);
             auto_next_update[lev] = -1.0;
             auto_ready[lev] = 0;
             // Regridding can leave vectors at the same size; resize alone would retain
             // clocks/amplitudes belonging to a different set of boxes.
             pb_interval[lev].assign(pb_ba[lev].size(), -one);
             pb_local_etime[lev].assign(pb_ba[lev].size(), zero);
             pb_nupdate[lev].assign(pb_ba[lev].size(), 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

Store the integer perturbation-type selector for one AMR level.

Parameters

lev	AMR level index.
pert_type	Perturbation type selected for this level.
max_level	Maximum AMR level configured for the run.

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

◆ netZeroBuoyant()

void TurbulentPerturbation::netZeroBuoyant ( const int & lev )

inline

Remove the mean of the box-perturbation field, so the net buoyant forcing is zero.

Each perturbation box's mean is summed over every grid box it overlaps and over all ranks, so the correction does not depend on the domain decomposition.

Parameters

lev	AMR level index.

     {
         const int nboxes = pb_ba[lev].size();
         if (nboxes == 0) { return; }
  
         // Sum and cell count of each perturbation box, then one reduction over the ranks
         amrex::Vector<amrex::Real> sums(2*nboxes, zero);
         for (amrex::MFIter mfi(pb_cell[lev]); mfi.isValid(); ++mfi) {
             const amrex::Box& vbx = mfi.validbox();
             const amrex::Array4<const amrex::Real>& pert_cell = pb_cell[lev].const_array(mfi);
             for (const auto& isect : pb_ba[lev].intersections(vbx)) {
                 sums[2*isect.first  ] += box_sum(pert_cell, isect.second);
                 sums[2*isect.first+1] += auto_mode
                     ? box_sum(auto_mask[lev].const_array(mfi), isect.second)
                     : static_cast<amrex::Real>(isect.second.numPts());
             }
         }
         amrex::ParallelAllReduce::Sum(sums.data(), static_cast<int>(sums.size()),
                                       amrex::ParallelContext::CommunicatorSub());
  
         // Normalizing the adjustment based on how many boxes there are; each box's
         // value is already normalized by the number of cells within it
         tpi_net_buoyant = zero;
         for (int boxIdx = 0; boxIdx < nboxes; boxIdx++) {
             pb_netZero[lev][boxIdx] = (sums[2*boxIdx+1] > zero) ? sums[2*boxIdx] / sums[2*boxIdx+1] : zero;
             tpi_net_buoyant += pb_netZero[lev][boxIdx];
         }
         tpi_pert_adjust = tpi_net_buoyant / static_cast<amrex::Real>(nboxes);
         if (auto_mode) {
             amrex::Real sum = zero, count = zero;
             for (int b = 0; b < nboxes; ++b) { sum += sums[2*b]; count += sums[2*b+1]; }
             tpi_pert_adjust = count > zero ? sum / count : zero;
         }
  
         const amrex::Real adjust = tpi_pert_adjust;
         for (amrex::MFIter mfi(pb_cell[lev]); mfi.isValid(); ++mfi) {
             const amrex::Box& vbx = mfi.validbox();
             const amrex::Array4<amrex::Real>& pert_cell = pb_cell[lev].array(mfi);
             const bool use_mask = auto_mode;
             const auto mask = use_mask ? auto_mask[lev].const_array(mfi)
                                        : amrex::Array4<const amrex::Real>{};
             for (const auto& isect : pb_ba[lev].intersections(vbx)) {
                 ParallelFor(isect.second, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
                     if (!use_mask || mask(i,j,k) > zero) { 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

Assign pseudo-random perturbations to cells in one perturbation box.

Parameters

boxIdx	Perturbation-box index on the level.
lev	AMR level index.
m_ixtype	Index type used to align perturbation storage.

     {
         for (amrex::MFIter mfi(pb_cell[lev],TileNoZ()); mfi.isValid(); ++mfi) {
             amrex::Box vbx = mfi.tilebox();
             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);
                 const bool use_mask = auto_mode;
                 const auto mask = use_mask ? auto_mask[lev].const_array(mfi)
                                            : amrex::Array4<const amrex::Real>{};
  
                 if (pt_type[lev] == 2  || pt_type[lev] == 3) { // CPM: one value per perturbation box
                     const amrex::Real rand_number_const = erf_hash_rng::hash_symmetric(
                         boxIdx, pb_nupdate[lev][boxIdx], 0, hash_comp_cpm, lev, tpi_seed);
                     ParallelFor(ubx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
                         pert_cell(i,j,k) = (use_mask && mask(i,j,k) == zero) ? zero : rand_number_const * amp_copy;
                     });
                 } else { // one value per cell, keyed by its global index and the box's refresh count
                     const std::uint64_t cell_seed = erf_hash_rng::splitmix64(
                         tpi_seed ^ erf_hash_rng::splitmix64(static_cast<std::uint64_t>(pb_nupdate[lev][boxIdx])));
                     ParallelFor(ubx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
                         pert_cell(i,j,k) = (use_mask && mask(i,j,k) == zero) ? zero
                             : erf_hash_rng::hash_symmetric(i, j, k, hash_comp_cell, lev, cell_seed) * amp_copy;
                     });
                 }
             }
         }
     }

Referenced by calc_tpi_uniform_update(), and calc_tpi_update().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ set_uniform_advection()

void TurbulentPerturbation::set_uniform_advection	(	int	lev,
		double	time,
		amrex::Real	speed,
		double	wind_period
	)

inline

     {
         AMREX_ALWAYS_ASSERT(uniform_refresh_enabled(lev));
         AMREX_ALWAYS_ASSERT(speed >= zero && std::isfinite(speed));
         auto& clock = uniform_clock[lev];
         clock.interval = speed > zero
             ? double(tpi_layers) * double(std::min(tpi_Lpb[lev], tpi_Wpb[lev])) / double(speed)
             : 0.0;
         clock.next_refresh = time;
         clock.next_wind_update = time + wind_period;
         clock.apply = false;
     }

Here is the call graph for this function:

◆ uniform_refresh_enabled()

bool TurbulentPerturbation::uniform_refresh_enabled ( int lev ) const

inline

     {
         return uniform_mode && (pt_type[lev] == 2 || pt_type[lev] == 3);
     }

Referenced by apply_tpi(), calc_tpi_uniform_update(), init_tpi(), and set_uniform_advection().

Here is the caller graph for this function:

◆ uniform_refresh_interval()

double TurbulentPerturbation::uniform_refresh_interval ( int lev ) const

inline

438 { return uniform_clock[lev].interval; }

◆ uniform_sampling_face()

bool TurbulentPerturbation::uniform_sampling_face ( int face ) const

inline

     {
         return tpi_direction[face < 2 ? face : face+1] != 0;
     }

◆ uniform_sampling_height()

amrex::Real TurbulentPerturbation::uniform_sampling_height ( ) const

inline

415 { return uniform_height; }

◆ uniform_sampling_period()

amrex::Real TurbulentPerturbation::uniform_sampling_period ( ) const

inline

416 { return uniform_wind_period; }

◆ uniform_wind_update_due()

bool TurbulentPerturbation::uniform_wind_update_due	(	int	lev,
		double	time
	)		const

inline

     {
         return time >= uniform_clock[lev].next_wind_update;
     }

◆ zero_amp()

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

inline

Reset CPM perturbation amplitudes in one perturbation box to zero.

Parameters

boxIdx	Perturbation-box index on the level.
lev	AMR level index.
m_ixtype	Index type used to align perturbation storage.

     {
  
         for (amrex::MFIter mfi(pb_cell[lev],TileNoZ()); mfi.isValid(); ++mfi) {
             amrex::Box vbx = mfi.tilebox();
             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

◆ auto_mask

amrex::Vector<amrex::MultiFab> TurbulentPerturbation::auto_mask

private

Referenced by configure_automatic(), init_tpi(), netZeroBuoyant(), and pseudoRandomPert().

◆ auto_mode

bool TurbulentPerturbation::auto_mode {false}

private

Referenced by calc_tpi_amp(), calc_tpi_uniform_update(), calc_tpi_update(), init_tpi(), netZeroBuoyant(), and pseudoRandomPert().

◆ auto_next_update

amrex::Vector<double> TurbulentPerturbation::auto_next_update

private

Referenced by automatic_update_due(), configure_automatic(), and init_tpi().

◆ auto_ready

amrex::Vector<int> TurbulentPerturbation::auto_ready

private

Referenced by automatic_update_due(), calc_tpi_uniform_update(), calc_tpi_update(), configure_automatic(), and init_tpi().

◆ hash_comp_cell

TurbulentPerturbation::hash_comp_cell = 0

staticconstexprprivate

Hash salt of the per-cell perturbation values.

Referenced by pseudoRandomPert().

◆ hash_comp_cpm

TurbulentPerturbation::hash_comp_cpm = 1

staticconstexprprivate

Hash salt of the per-box CPM values.

Referenced by pseudoRandomPert().

◆ hash_comp_interval

TurbulentPerturbation::hash_comp_interval = 2

staticconstexprprivate

Hash salt of the refresh-interval variation.

Referenced by calc_tpi_update().

◆ input_Ug

TurbulentPerturbation::input_Ug

private

Input geostrophic wind speed used to scale CPM perturbations.

Referenced by calc_tpi_amp(), and init_tpi().

◆ input_w_amp

TurbulentPerturbation::input_w_amp

private

Input vertical velocity perturbation amplitude.

Referenced by calc_tpi_amp(), and init_tpi().

◆ pb_amp

TurbulentPerturbation::pb_amp

private

Perturbation amplitude for each perturbation box.

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

◆ pb_ba

TurbulentPerturbation::pb_ba

Perturbation-box BoxArray for each AMR level.

Referenced by apply_tpi(), calc_tpi_meanMag(), calc_tpi_uniform_update(), calc_tpi_update(), configure_automatic(), init_tpi(), netZeroBuoyant(), pseudoRandomPert(), and zero_amp().

◆ pb_cell

TurbulentPerturbation::pb_cell

Per-cell perturbation amplitude storage.

Referenced by calc_tpi_uniform_update(), configure_automatic(), init_tpi(), make_sources(), netZeroBuoyant(), pseudoRandomPert(), and zero_amp().

◆ pb_dir

TurbulentPerturbation::pb_dir

Mean velocity direction for each perturbation box.

Referenced by calc_tpi_meanMag(), calc_tpi_update(), configure_automatic(), and init_tpi().

◆ pb_interval

TurbulentPerturbation::pb_interval

private

Perturbation update interval for each box [s].

Referenced by calc_tpi_uniform_update(), calc_tpi_update(), configure_automatic(), and init_tpi().

◆ pb_local_etime

TurbulentPerturbation::pb_local_etime

private

Local elapsed time for each perturbation box [s].

Referenced by calc_tpi_uniform_update(), calc_tpi_update(), configure_automatic(), and init_tpi().

◆ pb_mag

TurbulentPerturbation::pb_mag

Mean velocity magnitude for each perturbation box [m/s].

Referenced by calc_tpi_amp(), calc_tpi_meanMag(), calc_tpi_update(), configure_automatic(), and init_tpi().

◆ pb_netZero

TurbulentPerturbation::pb_netZero

private

Per-box storage used for net-zero buoyancy calculation.

Referenced by init_tpi(), and netZeroBuoyant().

◆ pb_nupdate

TurbulentPerturbation::pb_nupdate

private

Number of times each perturbation box has been refreshed; keys its random numbers.

Referenced by calc_tpi_uniform_update(), calc_tpi_update(), init_tpi(), and pseudoRandomPert().

◆ pb_Ug

amrex::Vector<amrex::Vector<amrex::Real> > TurbulentPerturbation::pb_Ug

private

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

◆ perturbation_khi

TurbulentPerturbation::perturbation_khi

private

Upper vertical index bound for perturbations.

Referenced by init_tpi().

◆ perturbation_klo

TurbulentPerturbation::perturbation_klo

private

Lower vertical index bound for perturbations.

Referenced by init_tpi().

◆ pt_type

TurbulentPerturbation::pt_type

Integer perturbation type for each AMR level.

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

◆ ref_ratio

TurbulentPerturbation::ref_ratio

private

Refinement ratios used by multilevel perturbation scaling.

Referenced by calc_tpi_amp(), and init_tpi().

◆ tpi_boxDim

TurbulentPerturbation::tpi_boxDim

private

Dimensions of each perturbation box.

Referenced by automatic_face_box(), and init_tpi().

◆ tpi_direction

TurbulentPerturbation::tpi_direction

private

Boundary directions where perturbations are applied.

Referenced by init_tpi(), and uniform_sampling_face().

◆ tpi_Hpb

TurbulentPerturbation::tpi_Hpb

private

Perturbation-box height for each AMR level [m].

Referenced by calc_tpi_amp(), and init_tpi().

◆ tpi_layers

TurbulentPerturbation::tpi_layers

private

Number of layers of perturbation boxes.

Referenced by automatic_face_box(), calc_tpi_uniform_update(), calc_tpi_update(), init_tpi(), and set_uniform_advection().

◆ tpi_Lpb

TurbulentPerturbation::tpi_Lpb

private

Perturbation-box length for each AMR level [m].

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

◆ tpi_lref

TurbulentPerturbation::tpi_lref

private

Perturbation-box reference length for each AMR level [m].

Referenced by calc_tpi_update(), and init_tpi().

◆ tpi_net_buoyant

TurbulentPerturbation::tpi_net_buoyant

private

Accumulated net buoyant perturbation used for correction.

Referenced by init_tpi(), and netZeroBuoyant().

◆ tpi_nonDim

TurbulentPerturbation::tpi_nonDim

private

Nondimensional number used by the perturbation formulation.

Referenced by calc_tpi_amp(), and init_tpi().

◆ tpi_offset

TurbulentPerturbation::tpi_offset

private

Cell offset for the start of the perturbation region.

Referenced by automatic_face_box(), and init_tpi().

◆ tpi_pert_adjust

TurbulentPerturbation::tpi_pert_adjust

private

Per-cell perturbation adjustment used for net-zero buoyancy.

Referenced by init_tpi(), and netZeroBuoyant().

◆ tpi_seed

TurbulentPerturbation::tpi_seed = 0

private

Seed of the hashed random numbers (1024 with erf.fix_random_seed, else erf.random_seed, else 0).

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

◆ tpi_Ti

TurbulentPerturbation::tpi_Ti

private

Temperature intensity used by the perturbation formulation.

Referenced by calc_tpi_amp(), and init_tpi().

◆ tpi_Tinf

TurbulentPerturbation::tpi_Tinf

private

Reference temperature used by the perturbation formulation [K].

Referenced by calc_tpi_amp(), and init_tpi().

◆ tpi_Wpb

TurbulentPerturbation::tpi_Wpb

private

Perturbation-box width for each AMR level [m].

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

◆ uniform_clock

amrex::Vector<UniformClock> TurbulentPerturbation::uniform_clock

private

Referenced by apply_tpi(), calc_tpi_uniform_update(), init_tpi(), set_uniform_advection(), uniform_refresh_interval(), and uniform_wind_update_due().

◆ uniform_height

amrex::Real TurbulentPerturbation::uniform_height {amrex::Real(20)}

private

Referenced by init_tpi(), and uniform_sampling_height().

◆ uniform_mode

bool TurbulentPerturbation::uniform_mode {true}

private

Referenced by init_tpi(), and uniform_refresh_enabled().

◆ uniform_wind_period

amrex::Real TurbulentPerturbation::uniform_wind_period {amrex::Real(600)}

private

Referenced by init_tpi(), and uniform_sampling_period().

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

Source/DataStructs/ERF_TurbPertStruct.H

Classes

Public Member Functions

Static Public Member Functions

Public Attributes

Private Attributes

Static Private Attributes

Detailed Description

Constructor & Destructor Documentation

◆ ~TurbulentPerturbation()

Member Function Documentation

◆ apply_tpi()

◆ automatic_face_box()

◆ automatic_update_due()

◆ box_sum()

◆ calc_tpi_amp()

◆ calc_tpi_meanMag()

◆ calc_tpi_uniform_update()

◆ calc_tpi_update()

◆ configure_automatic()

◆ debug()

◆ init_tpi()

◆ init_tpi_type()

◆ netZeroBuoyant()

◆ pseudoRandomPert()

◆ set_uniform_advection()

◆ uniform_refresh_enabled()

◆ uniform_refresh_interval()

◆ uniform_sampling_face()

◆ uniform_sampling_height()

◆ uniform_sampling_period()

◆ uniform_wind_update_due()

◆ zero_amp()

Member Data Documentation

◆ auto_mask

◆ auto_mode

◆ auto_next_update

◆ auto_ready

◆ hash_comp_cell

◆ hash_comp_cpm

◆ hash_comp_interval

◆ input_Ug

◆ input_w_amp

◆ pb_amp

◆ pb_ba

◆ pb_cell

◆ pb_dir

◆ pb_interval

◆ pb_local_etime

◆ pb_mag

◆ pb_netZero

◆ pb_nupdate

◆ pb_Ug

◆ perturbation_khi

◆ perturbation_klo

◆ pt_type

◆ ref_ratio

◆ tpi_boxDim

◆ tpi_direction

◆ tpi_Hpb

◆ tpi_layers

◆ tpi_Lpb

◆ tpi_lref

◆ tpi_net_buoyant

◆ tpi_nonDim

◆ tpi_offset

◆ tpi_pert_adjust

◆ tpi_seed

◆ tpi_Ti

◆ tpi_Tinf

◆ tpi_Wpb

◆ uniform_clock

◆ uniform_height

◆ uniform_mode

◆ uniform_wind_period