ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
TurbulentPerturbation Struct Reference

#include <ERF_TurbPertStruct.H>

Collaboration diagram for TurbulentPerturbation:

Public Member Functions

 ~TurbulentPerturbation ()
 
void init_tpi (const int lev, const amrex::IntVect &nx, const amrex::GpuArray< amrex::Real, 3 > dx, const amrex::BoxArray &ba, const amrex::DistributionMapping &dm, const int ngrow_state)
 
void calc_tpi_update (const int lev, const amrex::Real 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 &boxIdx, const amrex::Real &interval)
 
void pseudoRandomPert (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)
 
amrex::Real * get_pb_mag ()
 
amrex::Real * get_pb_netZero ()
 
void debug (amrex::Real)
 

Public Attributes

std::string pp_prefix {"erf"}
 
int pt_type
 
amrex::Vector< amrex::BoxArray > pb_ba
 
amrex::Vector< amrex::Real > pb_mag
 
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::Real tpi_Hpb
 
amrex::Real tpi_Lpb
 
amrex::Real tpi_Wpb
 
amrex::Real tpi_lref
 
amrex::Real tpi_net_buoyant
 
amrex::Real tpi_pert_adjust
 
amrex::Vector< amrex::Real > pb_interval
 
amrex::Vector< amrex::Real > pb_local_etime
 
amrex::Vector< amrex::Real > pb_amp
 
amrex::Vector< amrex::Real > pb_netZero
 

Detailed Description

Container holding quantities related to turbulent perturbation parameters

Constructor & Destructor Documentation

◆ ~TurbulentPerturbation()

TurbulentPerturbation::~TurbulentPerturbation ( )
inline
22 {}

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
251  {
252  for (int boxIdx = 0; boxIdx < pb_ba[lev].size(); boxIdx++) {
253  amrex::Box pbx = amrex::convert(pb_ba[lev][boxIdx], m_ixtype);
254  amrex::Box ubx = pbx & vbx;
255  if (ubx.ok()) {
256  ParallelFor(ubx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
257  src_arr(i,j,k,comp) += pert_cell(i,j,k);
258 
259  // For box region debug only
260  #ifdef INDEX_PERTURB
261  src_arr(i,j,k,comp) = (amrex::Real) (boxIdx + 5.);
262  #endif
263  });
264  }
265  }
266  }
amrex::Vector< amrex::BoxArray > pb_ba
Definition: ERF_TurbPertStruct.H:531

Referenced by make_sources().

Here is the caller graph for this function:

◆ calc_tpi_amp()

void TurbulentPerturbation::calc_tpi_amp ( const int &  boxIdx,
const amrex::Real &  interval 
)
inline
271  {
272  amrex::Real Um = pb_mag[boxIdx];
273  pb_amp[boxIdx] = 0.; // Safety step
274 
275  amrex::Real beta = 1./tpi_Tinf; // Thermal expansion coefficient
276 
277  // Pseudo Random temperature the ignores scale when mechanically tripping turbulence
278  amrex::Real g = CONST_GRAV;
279  if (tpi_Ti > 0.) g = (tpi_nonDim * Um * Um) / (tpi_Ti * tpi_Hpb);
280 
281  // Ma and Senocak (2023) Eq. 8, solving for delta phi
282  pb_amp[boxIdx] = (tpi_nonDim * Um * Um) / (g * beta * tpi_Hpb);
283 
284  if (pt_type == 0) {
285  // Performing this step converts the perturbation proportionality into
286  // the forcing term
287  // Ma & Senocak (2023) Eq. 7
288  pb_amp[boxIdx] /= interval;
289  }
290  }
constexpr amrex::Real CONST_GRAV
Definition: ERF_Constants.H:21
amrex::Real tpi_Tinf
Definition: ERF_TurbPertStruct.H:550
amrex::Vector< amrex::Real > pb_amp
Definition: ERF_TurbPertStruct.H:564
amrex::Real tpi_Ti
Definition: ERF_TurbPertStruct.H:549
amrex::Vector< amrex::Real > pb_mag
Definition: ERF_TurbPertStruct.H:532
amrex::Real tpi_nonDim
Definition: ERF_TurbPertStruct.H:548
int pt_type
Definition: ERF_TurbPertStruct.H:528
amrex::Real tpi_Hpb
Definition: ERF_TurbPertStruct.H:553

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
386  {
387  // Creating local copy of PB box array and magnitude
388  const amrex::BoxArray m_pb_ba = pb_ba[lev];
389  amrex::Real* m_pb_mag = get_pb_mag();
390  m_pb_mag[boxIdx] = 0.; // Safety step
391 
392  // Storage of averages per PB
393  // Index: 0=u (vol/slab_lo), 1=v (vol/slab_lo)
394  // 2=u (slab_hi), 3=v (slab_hi)
395  int n_avg = 4;
396  amrex::Vector<amrex::Real> avg_h(n_avg,0.);
397  amrex::Gpu::DeviceVector<amrex::Real> avg_d(n_avg,0.);
398  amrex::Real* avg = avg_d.data();
399 
400  // Averaging u & v components in single MFIter
401  for (amrex::MFIter mfi(mf_cons, TileNoZ()); mfi.isValid(); ++mfi) {
402 
403  // CC valid box (inherited from mf_cons)
404  const amrex::Box& vbx = mfi.validbox();
405 
406  // Box logic for u velocity
407  auto ixtype_u = mf_xvel.boxArray().ixType();
408  amrex::Box vbx_u = amrex::convert(vbx,ixtype_u);
409  amrex::Box pbx_u = amrex::convert(m_pb_ba[boxIdx], ixtype_u);
410  amrex::Box ubx_u = pbx_u & vbx_u;
411 
412  // Box logic for v velocity
413  auto ixtype_v = mf_yvel.boxArray().ixType();
414  amrex::Box vbx_v = amrex::convert(vbx,ixtype_v);
415  amrex::Box pbx_v = amrex::convert(m_pb_ba[boxIdx], ixtype_v);
416  amrex::Box ubx_v = pbx_v & vbx_v;
417 
418  // Operation over box union (U)
419  if (ubx_u.ok()) {
420  const amrex::Array4<const amrex::Real>& xvel_arry = mf_xvel.const_array(mfi);
421 
422  #ifdef USE_VOLUME_AVERAGE
423  int npts = ubx_u.numPts();
424  amrex::Real norm = 1.0 / (amrex::Real) npts;
425  ParallelFor(amrex::Gpu::KernelInfo().setReduction(true), ubx_u, [=]
426  AMREX_GPU_DEVICE(int i, int j, int k, amrex::Gpu::Handler const& handler) noexcept {
427  amrex::Gpu::deviceReduceSum(&avg[0], xvel_arry(i,j,k)*norm, handler);
428  });
429  #endif // USE_VOLUME_AVERAGE
430 
431  #ifdef USE_SLAB_AVERAGE
432  amrex::Box ubxSlab_lo = makeSlab(ubx_u,2,ubx_u.smallEnd(2));
433  amrex::Box ubxSlab_hi = makeSlab(ubx_u,2,ubx_u.bigEnd(2));
434  int npts_lo = ubxSlab_lo.numPts();
435  int npts_hi = ubxSlab_hi.numPts();
436  amrex::Real norm_lo = 1.0 / (amrex::Real) npts_lo;
437  amrex::Real norm_hi = 1.0 / (amrex::Real) npts_hi;
438 
439  // Average u in the low slab
440  ParallelFor(amrex::Gpu::KernelInfo().setReduction(true), ubxSlab_lo, [=]
441  AMREX_GPU_DEVICE(int i, int j, int k, amrex::Gpu::Handler const& handler) noexcept {
442  amrex::Gpu::deviceReduceSum(&avg[0], xvel_arry(i,j,k)*norm_lo, handler);
443  });
444 
445  // Average u in the high slab
446  ParallelFor(amrex::Gpu::KernelInfo().setReduction(true), ubxSlab_hi, [=]
447  AMREX_GPU_DEVICE(int i, int j, int k, amrex::Gpu::Handler const& handler) noexcept {
448  amrex::Gpu::deviceReduceSum(&avg[2], xvel_arry(i,j,k)*norm_hi, handler);
449  });
450  #endif // USE_SLAB_AVERAGE
451  } // if
452 
453  // Operation over box union (V)
454  if (ubx_v.ok()) {
455  const amrex::Array4<const amrex::Real>& yvel_arry = mf_yvel.const_array(mfi);
456 
457  #ifdef USE_VOLUME_AVERAGE
458  int npts = ubx_v.numPts();
459  amrex::Real norm = 1.0 / (amrex::Real) npts;
460  ParallelFor(amrex::Gpu::KernelInfo().setReduction(true), ubx_v, [=]
461  AMREX_GPU_DEVICE(int i, int j, int k, amrex::Gpu::Handler const& handler) noexcept {
462  amrex::Gpu::deviceReduceSum(&avg[1], yvel_arry(i,j,k)*norm, handler);
463  });
464  #endif // USE_VOLUME_AVERAGE
465 
466  #ifdef USE_SLAB_AVERAGE
467  amrex::Box ubxSlab_lo = makeSlab(ubx_v,2,ubx_v.smallEnd(2));
468  amrex::Box ubxSlab_hi = makeSlab(ubx_v,2,ubx_v.bigEnd(2));
469  int npts_lo = ubxSlab_lo.numPts();
470  int npts_hi = ubxSlab_hi.numPts();
471  amrex::Real norm_lo = 1.0 / (amrex::Real) npts_lo;
472  amrex::Real norm_hi = 1.0 / (amrex::Real) npts_hi;
473 
474  // Average v in the low slab
475  ParallelFor(amrex::Gpu::KernelInfo().setReduction(true), ubxSlab_lo, [=]
476  AMREX_GPU_DEVICE(int i, int j, int k, amrex::Gpu::Handler const& handler) noexcept {
477  amrex::Gpu::deviceReduceSum(&avg[1], yvel_arry(i,j,k)*norm_lo, handler);
478  });
479 
480  // Average v in the high slab
481  ParallelFor(amrex::Gpu::KernelInfo().setReduction(true), ubxSlab_hi, [=]
482  AMREX_GPU_DEVICE(int i, int j, int k, amrex::Gpu::Handler const& handler) noexcept {
483  amrex::Gpu::deviceReduceSum(&avg[3], yvel_arry(i,j,k)*norm_hi, handler);
484  });
485  #endif // USE_SLAB_AVERAGE
486  } // if
487  } // MFIter
488 
489  // Copy from device back to host
490  amrex::Gpu::copy(amrex::Gpu::deviceToHost, avg_d.begin(), avg_d.end(), avg_h.begin());
491 
492  // Computing the average magnitude within PB
493  #ifdef USE_VOLUME_AVERAGE
494  m_pb_mag[boxIdx] = sqrt(avg_h[0]*avg_h[0] + avg_h[1]*avg_h[1]);
495  #endif
496 
497  #ifdef USE_SLAB_AVERAGE
498  m_pb_mag[boxIdx] = 0.5*( sqrt(avg_h[0]*avg_h[0] + avg_h[1]*avg_h[1])
499  + sqrt(avg_h[2]*avg_h[2] + avg_h[3]*avg_h[3]));
500  #endif
501  }
AMREX_FORCE_INLINE amrex::IntVect TileNoZ()
Definition: ERF_TileNoZ.H:11
amrex::Real * get_pb_mag()
Definition: ERF_TurbPertStruct.H:504

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 amrex::Real  dt,
amrex::MultiFab &  mf_xvel,
amrex::MultiFab &  mf_yvel,
amrex::MultiFab &  mf_cons 
)
inline
189  {
190  // Resettubg the net buoyant force value
191  tpi_net_buoyant = 0.;
192 
193  // Setting random number generator for update interval
194  srand( (unsigned) time(NULL) );
195 
196  auto m_ixtype = mf_cons.boxArray().ixType(); // safety step
197 
198  for (int boxIdx = 0; boxIdx < pb_ba[lev].size(); boxIdx++) {
199 
200  // Check if the local elapsed time is greater than the update interval
201  if ( pb_interval[boxIdx] <= pb_local_etime[boxIdx] ) {
202 
203  // Compute mean velocity of each perturbation box
204  calc_tpi_meanMag_perBox(boxIdx, lev, mf_cons, mf_xvel, mf_yvel);
205 
206  // Only the rank owning the box will be able to access the storage location
207  // Done for parallelism to avoid Inf being stored in array
208  if (pb_mag[boxIdx] !=0.) {
209  amrex::Real interval = tpi_lref / pb_mag[boxIdx];
210  pb_interval[boxIdx] = RandomReal(0.9*interval,1.1*interval); // 10% variation
211  }
212 
213  // Trigger amplitude calculation per perturbation box
214  calc_tpi_amp(boxIdx, pb_interval[boxIdx]);
215 
216  // Trigger random amplitude storage per cell within perturbation box
217  pseudoRandomPert(boxIdx, lev, m_ixtype);
218 
219  // Reset local elapsed time
220  pb_local_etime[boxIdx] = 0.;
221  } else {
222  // Increase by timestep of level 0
223  pb_local_etime[boxIdx] += dt;
224  } // if
225 
226  // Per iteration operation of net-zero buoyant force check
227  if (pb_mag[boxIdx] !=0.) netZeroBuoyantAdd(boxIdx, lev);
228  tpi_net_buoyant += pb_netZero[boxIdx];
229  } // for
230 
231  // Normalizing the adjustment based on how many boxes there are
232  // the values within the array is already normalized by the number
233  // of cells within each box
234  tpi_pert_adjust = tpi_net_buoyant / (amrex::Real) pb_ba[lev].size();
235 
236  // Per iteration operation of net-zero buoyant force adjustment
237  for (int boxIdx = 0; boxIdx < pb_ba[lev].size(); boxIdx++) {
238  if (pb_mag[boxIdx] !=0.) netZeroBuoyantAdjust(boxIdx, lev);
239  }
240  }
void pseudoRandomPert(const int &boxIdx, const int &lev, const amrex::IndexType &m_ixtype)
Definition: ERF_TurbPertStruct.H:295
amrex::Vector< amrex::Real > pb_local_etime
Definition: ERF_TurbPertStruct.H:563
void netZeroBuoyantAdjust(const int &boxIdx, const int &lev)
Definition: ERF_TurbPertStruct.H:357
void calc_tpi_amp(const int &boxIdx, const amrex::Real &interval)
Definition: ERF_TurbPertStruct.H:269
amrex::Vector< amrex::Real > pb_netZero
Definition: ERF_TurbPertStruct.H:565
amrex::Real tpi_pert_adjust
Definition: ERF_TurbPertStruct.H:559
amrex::Vector< amrex::Real > pb_interval
Definition: ERF_TurbPertStruct.H:562
amrex::Real tpi_lref
Definition: ERF_TurbPertStruct.H:556
void calc_tpi_meanMag_perBox(const int &boxIdx, const int &lev, amrex::MultiFab &mf_cons, amrex::MultiFab &mf_xvel, amrex::MultiFab &mf_yvel)
Definition: ERF_TurbPertStruct.H:380
void netZeroBuoyantAdd(const int &boxIdx, const int &lev)
Definition: ERF_TurbPertStruct.H:322
amrex::Real RandomReal(const amrex::Real min, const amrex::Real max)
Definition: ERF_TurbPertStruct.H:568
amrex::Real tpi_net_buoyant
Definition: ERF_TurbPertStruct.H:558
Here is the call graph for this function:

◆ debug()

void TurbulentPerturbation::debug ( amrex::Real  )
inline
509  {
510  /*
511  amrex::PrintToFile("BoxPerturbationOutput") << "#################### PB output at time = "
512  << time << " ####################\n";
513  amrex::PrintToFile("BoxPerturbationOutput") << " Using type: " << pt_type << "\n";
514  amrex::PrintToFile("BoxPerturbationOutput") << " Net: " << tpi_net_buoyant << " Adjust : " << tpi_pert_adjust << "\n";
515  for (int i = 0; i < pb_mag.size(); i++) {
516  amrex::PrintToFile("BoxPerturbationOutput") << "[" << i
517  << "] pb_Umag=" << pb_mag[i]
518  << " | pb_interval=" << pb_interval[i]
519  << " (" << pb_local_etime[i]
520  << ") | pb_amp=" << pb_amp[i] << "\n";
521  }
522  amrex::PrintToFile("BoxPerturbationOutput") << "\n";
523  */
524  }

◆ get_pb_mag()

amrex::Real* TurbulentPerturbation::get_pb_mag ( )
inline
504 { return pb_mag.data(); }

Referenced by calc_tpi_meanMag_perBox().

Here is the caller graph for this function:

◆ get_pb_netZero()

amrex::Real* TurbulentPerturbation::get_pb_netZero ( )
inline
505 { return pb_netZero.data(); }

Referenced by netZeroBuoyantAdd().

Here is the caller graph for this function:

◆ init_tpi()

void TurbulentPerturbation::init_tpi ( const int  lev,
const amrex::IntVect &  nx,
const amrex::GpuArray< amrex::Real, 3 >  dx,
const amrex::BoxArray &  ba,
const amrex::DistributionMapping &  dm,
const int  ngrow_state 
)
inline
36  {
37  // Initialization for some 0 dependent terms
38  tpi_Ti = 0.;
39  tpi_Tinf = 300.;
40 
41  amrex::ParmParse pp(pp_prefix);
42 
43  // Reading inputs, and placing assertion for the perturbation inflow to work
44  pp.queryarr("perturbation_box_dims",tpi_boxDim);
45  pp.queryarr("perturbation_direction",tpi_direction);
46  pp.query("perturbation_layers",tpi_layers);
47  pp.query("perturbation_offset",tpi_offset);
48 
49  pp.query("perturbation_nondimensional",tpi_nonDim);
50  pp.query("perturbation_T_infinity",tpi_Tinf);
51  pp.query("perturbation_T_intensity",tpi_Ti);
52 
53  // Check variables message
54  if (tpi_layers < 0) { amrex::Abort("Please provide a valid perturbation layer value (ie. 3-5)"); }
55  if (tpi_offset < 0) { amrex::Abort("Please provide a valid inflow cell offset value for perturbation region (ie. 0-5)"); }
56  for (int i = 0; i < tpi_boxDim.size(); i++) {
57  if (tpi_boxDim[i] == 0) { amrex::Abort("Please provide valid dimensions for perturbation boxes."); }
58  }
59  if (tpi_nonDim < 0.) { amrex::Abort("Please provide a valid nondimensional number (ie. Ri = 0.042)"); }
60  if (tpi_Tinf < 0.) { amrex::Abort("Please provide a valid ambient temperature value (ie. T_0 = T_infty)"); }
61  if (tpi_Ti < 0.) { amrex::Abort("Please provide a valid temperature intensity value (ie. 0-1.0)"); }
62 
63  // Creating perturbation regions and initializing with generic size. Temporary size for now
64  amrex::Box lo_x_bx(amrex::IntVect(0), amrex::IntVect(1), amrex::IntVect(0));
65  amrex::Box hi_x_bx(amrex::IntVect(0), amrex::IntVect(1), amrex::IntVect(0));
66  amrex::Box lo_y_bx(amrex::IntVect(0), amrex::IntVect(1), amrex::IntVect(0));
67  amrex::Box hi_y_bx(amrex::IntVect(0), amrex::IntVect(1), amrex::IntVect(0));
68 
69  // Create a temporary box list to accumulate all the perturbation regions after box modification
70  amrex::BoxList tmp_bl;
71 
72  // boxSize for individual boxes
73  amrex::IntVect boxSize(tpi_boxDim[0],tpi_boxDim[1],tpi_boxDim[2]);
74 
75  // Starting logic to set the size of the perturbation region(s)
76  //amrex::PrintToFile("BoxPerturbationOutput") << "Setting perturbation region in:";
77  // ***** X-direction perturbation *****
78  if (tpi_direction[0]) { // West
79  lo_x_bx.setSmall(amrex::IntVect(tpi_offset, tpi_direction[1]*tpi_offset, 0));
80  lo_x_bx.setBig(amrex::IntVect((tpi_layers*tpi_boxDim[0]-1)+tpi_offset, nx[1]-(tpi_direction[4]*tpi_offset), nx[2]));
81  //amrex::PrintToFile("BoxPerturbationOutput") << " West face";
82  }
83 
84  if (tpi_direction[3]) { // East
85  hi_x_bx.setSmall(amrex::IntVect(nx[0]-((tpi_layers*tpi_boxDim[0]-1)+tpi_offset), tpi_direction[1]*tpi_offset, 0));
86  hi_x_bx.setBig(amrex::IntVect(nx[0]-tpi_offset, nx[1]-(tpi_direction[4]*tpi_offset), nx[2]));
87  //amrex::PrintToFile("BoxPerturbationOutput") << " East face";
88  }
89 
90  // ***** Y-direction Perturbation *****
91  if (tpi_direction[1]) { // North
92  lo_y_bx.setSmall(amrex::IntVect(tpi_direction[0]*tpi_offset, tpi_offset, 0));
93  lo_y_bx.setBig(amrex::IntVect(nx[0]-tpi_direction[3]*tpi_offset, ((tpi_layers*tpi_boxDim[1])-1)+tpi_offset, nx[2]));
94  //amrex::PrintToFile("BoxPerturbationOutput") << " North face";
95  }
96 
97  if (tpi_direction[4]) { // South
98  hi_y_bx.setSmall(amrex::IntVect(tpi_direction[0]*tpi_offset, nx[1]-((tpi_layers*tpi_boxDim[1]-1)+tpi_offset), 0));
99  hi_y_bx.setBig(amrex::IntVect(nx[0]-tpi_direction[3]*tpi_offset, nx[1]-tpi_offset, nx[2]));
100  //amrex::PrintToFile("BoxPerturbationOutput") << " South face";
101  }
102 
103  if (tpi_direction[2] || tpi_direction[5]) { amrex::Abort("Currently not supporting z-direction flow perturbation"); }
104 
105  // Performing box union for intersecting perturbation regions to avoid overlapping sections (double counting at corners)
106  if (tpi_direction[0] && tpi_direction[1]) { // Reshaping South smallEnd
107  amrex::Box lo_x_lo_y_u = lo_x_bx & lo_y_bx;
108  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)));
109  }
110 
111  if (tpi_direction[3] && tpi_direction[1]) { // Reshaping South bigEnd
112  amrex::Box hi_x_lo_y_u = hi_x_bx & lo_y_bx;
113  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)));
114  }
115 
116  if (tpi_direction[0] && tpi_direction[4]) { // Reshaping North smallEnd
117  amrex::Box lo_x_hi_y_u = lo_x_bx & hi_y_bx;
118  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)));
119  }
120 
121  if (tpi_direction[3] && tpi_direction[4]) { // Reshaping North bigEnd
122  amrex::Box hi_x_hi_y_u = hi_x_bx & hi_y_bx;
123  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)));
124  }
125 
126  // Creating structure box array for conserved quantity
127  if (tpi_direction[0]) { tmp_bl.push_back(lo_x_bx); }
128  if (tpi_direction[1]) { tmp_bl.push_back(lo_y_bx); }
129  if (tpi_direction[3]) { tmp_bl.push_back(hi_x_bx); }
130  if (tpi_direction[4]) { tmp_bl.push_back(hi_y_bx); }
131  //amrex::PrintToFile("BoxPerturbationOutput") << "\nBoxList: " << tmp_bl << "\n";
132  amrex::BoxArray tmp_ba(tmp_bl);
133  tmp_ba.maxSize(boxSize);
134  pb_ba.push_back(tmp_ba);
135 
136  // Initializing mean magnitude vector
137  pb_mag.resize(pb_ba[lev].size(), 0.);
138  pb_netZero.resize(pb_ba[lev].size(), 0.);
139 
140  // Set size of vector and initialize
141  pb_interval.resize(pb_ba[lev].size(), -1.0);
142  pb_local_etime.resize(pb_ba[lev].size(), 0.0);
143  pb_amp.resize(pb_ba[lev].size(), 0.0);
144 
145  // Creating data array for perturbation amplitude storage
146  pb_cell.define(ba, dm, 1, ngrow_state);
147  pb_cell.setVal(0.);
148 
149  // Computing perturbation reference length
150  tpi_Lpb = tpi_boxDim[0]*dx[0];
151  tpi_Wpb = tpi_boxDim[1]*dx[1];
152  tpi_Hpb = tpi_boxDim[2]*dx[2];
154 
155  tpi_pert_adjust = 0.;
156  tpi_net_buoyant = 0.;
157 
158  /*
159  // Function check point message
160  amrex::PrintToFile("BoxPerturbationOutput") << "perturbation_box_dims: "
161  << tpi_boxDim[0] << " "
162  << tpi_boxDim[1] << " "
163  << tpi_boxDim[2] << "\n";
164  amrex::PrintToFile("BoxPerturbationOutput") << "perturbation_direction: "
165  << tpi_direction[0] << " "
166  << tpi_direction[1] << " "
167  << tpi_direction[2] << "\n\n";
168  amrex::PrintToFile("BoxPerturbationOutput") << "perturbation_layers: " << tpi_layers << "\n";
169  amrex::PrintToFile("BoxPerturbationOutput") << "perturbation_offset: " << tpi_offset << "\n\n";
170  amrex::PrintToFile("BoxPerturbationOutput") << "perturbation_nondimensional: " << tpi_nonDim << "\n";
171  amrex::PrintToFile("BoxPerturbationOutput") << "perturbation_T_infinity: " << tpi_Tinf << "\n";
172  amrex::PrintToFile("BoxPerturbationOutput") << "perturbation_T_intensity: " << tpi_Ti << "\n";
173  amrex::PrintToFile("BoxPerturbationOutput") << "Reference length per box = " << tpi_lref << "\n\n";
174  amrex::PrintToFile("BoxPerturbationOutput") << "Turbulent perturbation BoxArray:\n" << pb_ba[lev] << "\n";
175  */
176  }
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real pp(amrex::Real y)
Definition: ERF_MicrophysicsUtils.H:219
int tpi_layers
Definition: ERF_TurbPertStruct.H:541
amrex::Real tpi_Wpb
Definition: ERF_TurbPertStruct.H:555
amrex::Real tpi_Lpb
Definition: ERF_TurbPertStruct.H:554
std::string pp_prefix
Definition: ERF_TurbPertStruct.H:526
amrex::MultiFab pb_cell
Definition: ERF_TurbPertStruct.H:536
int tpi_offset
Definition: ERF_TurbPertStruct.H:542
amrex::Vector< int > tpi_boxDim
Definition: ERF_TurbPertStruct.H:544
amrex::Vector< int > tpi_direction
Definition: ERF_TurbPertStruct.H:545
Here is the call graph for this function:

◆ netZeroBuoyantAdd()

void TurbulentPerturbation::netZeroBuoyantAdd ( const int &  boxIdx,
const int &  lev 
)
inline
324  {
325  // Creating local copy of PB box array and magnitude
326  const amrex::BoxArray m_pb_ba = pb_ba[lev];
327  amrex::Real* m_pb_netZero = get_pb_netZero();
328 
329  // Create device array for summation
330  amrex::Vector<amrex::Real> avg_h(1,0.);
331  amrex::Gpu::DeviceVector<amrex::Real> avg_d(1,0.);
332  amrex::Real* avg = avg_d.data();
333 
334  // Iterates through the cells of each box and sum the white noise perturbation
335  for (amrex::MFIter mfi(pb_cell, TileNoZ()) ; mfi.isValid(); ++mfi) {
336  const amrex::Box& vbx = mfi.validbox();
337  amrex::Box pbx = amrex::convert(m_pb_ba[boxIdx], vbx.ixType());
338  amrex::Box ubx = pbx & vbx;
339  if (ubx.ok()) {
340  const amrex::Array4<const amrex::Real>& pert_cell = pb_cell.const_array(mfi);
341  int npts = ubx.numPts();
342  amrex::Real norm = 1.0 / (amrex::Real) npts;
343  ParallelFor(amrex::Gpu::KernelInfo().setReduction(true), ubx, [=]
344  AMREX_GPU_DEVICE(int i, int j, int k, amrex::Gpu::Handler const& handler) noexcept {
345  amrex::Gpu::deviceReduceSum(&avg[0], pert_cell(i,j,k)*norm, handler);
346  });
347  amrex::Gpu::copy(amrex::Gpu::deviceToHost, avg_d.begin(), avg_d.end(), avg_h.begin());
348 
349  // Assigning onto storage array
350  m_pb_netZero[boxIdx] = avg_h[0];
351  }
352  }
353  }
amrex::Real * get_pb_netZero()
Definition: ERF_TurbPertStruct.H:505

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
359  {
360  // Creating local copy of PB box array and magnitude
361  const amrex::BoxArray m_pb_ba = pb_ba[lev];
362  for (amrex::MFIter mfi(pb_cell, TileNoZ()) ; mfi.isValid(); ++mfi) {
363  const amrex::Box& vbx = mfi.validbox();
364  amrex::Box pbx = amrex::convert(m_pb_ba[boxIdx], vbx.ixType());
365  amrex::Box ubx = pbx & vbx;
366  if (ubx.ok()) {
367  const amrex::Real adjust = tpi_pert_adjust;
368  const amrex::Array4<amrex::Real>& pert_cell = pb_cell.array(mfi);
369  ParallelFor(ubx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
370  pert_cell(i,j,k) -= adjust;
371  });
372  }
373  }
374  }

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
298  {
299  // Seed the random generator at 1024UL for regression testing
300  int fix_random_seed = 0;
301  amrex::ParmParse pp("erf"); pp.query("fix_random_seed", fix_random_seed);
302  if (fix_random_seed) {
303  amrex::InitRandom(1024UL);
304  }
305 
306  for (amrex::MFIter mfi(pb_cell,TileNoZ()); mfi.isValid(); ++mfi) {
307  amrex::Box vbx = mfi.validbox();
308  amrex::Box pbx = amrex::convert(pb_ba[lev][boxIdx], m_ixtype);
309  amrex::Box ubx = pbx & vbx;
310  if (ubx.ok()) {
311  amrex::Real amp_copy = pb_amp[boxIdx];
312  const amrex::Array4<amrex::Real>& pert_cell = pb_cell.array(mfi);
313  ParallelForRNG(ubx, [=] AMREX_GPU_DEVICE(int i, int j, int k, const amrex::RandomEngine& engine) noexcept {
314  amrex::Real rand_double = amrex::Random(engine);
315  pert_cell(i,j,k) = (rand_double*2.0 - 1.0) * amp_copy;
316  });
317  }
318  }
319  }

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
569  {
570  amrex::Real r = (amrex::Real) rand() / (amrex::Real) RAND_MAX;
571  return min + r * (max - min);
572  }

Referenced by calc_tpi_update().

Here is the caller graph for this function:

Member Data Documentation

◆ pb_amp

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

◆ pb_ba

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

◆ pb_cell

amrex::MultiFab TurbulentPerturbation::pb_cell

◆ pb_interval

amrex::Vector<amrex::Real> TurbulentPerturbation::pb_interval
private

Referenced by calc_tpi_update(), and init_tpi().

◆ pb_local_etime

amrex::Vector<amrex::Real> TurbulentPerturbation::pb_local_etime
private

Referenced by calc_tpi_update(), and init_tpi().

◆ pb_mag

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

◆ pb_netZero

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

◆ pp_prefix

std::string TurbulentPerturbation::pp_prefix {"erf"}

Referenced by init_tpi().

◆ pt_type

int TurbulentPerturbation::pt_type

Referenced by calc_tpi_amp().

◆ 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::Real TurbulentPerturbation::tpi_Hpb
private

Referenced by calc_tpi_amp(), and init_tpi().

◆ tpi_layers

int TurbulentPerturbation::tpi_layers
private

Referenced by init_tpi().

◆ tpi_Lpb

amrex::Real TurbulentPerturbation::tpi_Lpb
private

Referenced by init_tpi().

◆ tpi_lref

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

◆ 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::Real TurbulentPerturbation::tpi_Wpb
private

Referenced by init_tpi().


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