ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
ERF_TurbPertStruct.H
Go to the documentation of this file.
1 #ifndef ERF_TURB_PERT_STRUCT_H_
2 #define ERF_TURB_PERT_STRUCT_H_
3 
4 #include <ERF_DataStruct.H>
5 #include <AMReX_MultiFabUtil.H>
6 #include <ERF_TileNoZ.H>
7 #include <time.h>
8 /**
9  * Container holding quantities related to turbulent perturbation parameters
10  */
11 
12 /* The general rule of thumb is to create a perturbation box size of 1/8th of the boundary layer height.
13  The boundary layer height can't be the height of the domain. The length and width of the box should
14  be twice the height of the box. If meandering flow is present, the width of the box should be take
15  the angle of the inflow into consideration.
16 */
17 
19 
20  public:
21 
23 
24  // Initializing Perturbation Region
25  // Currently only support perturbation in the x and y direction
26  // This portion of the code is only called once for initialization,
27  // therefore efficiency is not required. Dont spend too much time
28  // here other than to correctly setup the perturbation box region
29  void init_tpi (const int lev,
30  const amrex::IntVect& nx,
31  const amrex::GpuArray<amrex::Real,3> dx,
32  const amrex::BoxArray& ba,
33  const amrex::DistributionMapping& dm,
34  const int ngrow_state)
35 
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  }
177 
178 
179  // Perturbation update frequency check.
180  // This function will trigger calc_tpi_meanMag() and calc_tpi_amp(), when
181  // the local elapsed time is greater than the perturbation frequency.
182  // At each timestep a total buoyant force sanity check is made to ensure
183  // that the sum of the buoyant force introduced into the system is net-zero
184  void calc_tpi_update (const int lev,
185  const amrex::Real dt,
186  amrex::MultiFab& mf_xvel,
187  amrex::MultiFab& mf_yvel,
188  amrex::MultiFab& mf_cons)
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  }
241 
242  // Applying perturbation amplitude onto source term (Umphrey and Senocak 2016)
243  // This function does per cell application within the box union. Random perturbation
244  // is assigned in calc_tpi_update.
245  void apply_tpi (const int& lev,
246  const amrex::Box& vbx, // box union from upper level
247  const int& comp, // Component to modify
248  const amrex::IndexType& m_ixtype, // IntVect type of src_arr
249  const amrex::Array4<amrex::Real>& src_arr, // Array to apply perturbation
250  const amrex::Array4<amrex::Real const>& pert_cell)
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  }
267 
268  // Perturbation amplitude calculation
269  void calc_tpi_amp (const int& boxIdx,
270  const amrex::Real& interval)
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  }
291 
292  // Assigns pseudo-random (ie. white noise) perturbation to a storage cell, this
293  // value is then held constant for the duration of the update interval and assigned onto
294  // the source term
295  void pseudoRandomPert (const int& boxIdx,
296  const int& lev,
297  const amrex::IndexType& m_ixtype)
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  }
320 
321  // Checks for net-zero buoyant force introduction into the system
322  void netZeroBuoyantAdd (const int& boxIdx,
323  const int& lev)
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  }
354 
355  // If it's not a net-zero buoyant force, then adjust all cell by a normalized value
356  // to achieve this logic
357  void netZeroBuoyantAdjust (const int& boxIdx,
358  const int& lev)
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  }
375 
376 // TODO: Test the difference between these two for Source term perturbation
377 #define USE_VOLUME_AVERAGE
378  // Perturbation box mean velocity magnitude calculation
379  // This is pulled into the structure to also utilize during runtime
380  void calc_tpi_meanMag_perBox (const int& boxIdx,
381  const int& lev,
382  amrex::MultiFab& mf_cons,
383  amrex::MultiFab& mf_xvel,
384  amrex::MultiFab& mf_yvel)
385 
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  }
502 
503  // Quality of life function definitions
504  amrex::Real* get_pb_mag() { return pb_mag.data(); }
505  amrex::Real* get_pb_netZero() { return pb_netZero.data(); }
506 
507  // Output debug message into a file
508  void debug (amrex::Real /*time*/)
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  }
525 
526  std::string pp_prefix {"erf"};
527 
528  int pt_type; // TODO: May need to pass in solverChoice to replace this
529 
530  // Public data members
531  amrex::Vector<amrex::BoxArray> pb_ba; // PB box array
532  amrex::Vector<amrex::Real> pb_mag; // BP mean magnitude [m/s]
533 
534  // Perturbation amplitude cell storage
535  // This is after random assignment of equation (10) in Ma and Senocak 2023
536  amrex::MultiFab pb_cell;
537 
538  private:
539 
540  // Private data members
541  int tpi_layers; // Number of layers of perturbation boxes
542  int tpi_offset; // Cells to offset the start of the perturbation region
543 
544  amrex::Vector<int> tpi_boxDim; // Dimensions of each perturbation box
545  amrex::Vector<int> tpi_direction; // Direction of the peturbation
546 
547  // Richardson Formulation
548  amrex::Real tpi_nonDim; // Richardson number
549  amrex::Real tpi_Ti; // Temperature intensity value
550  amrex::Real tpi_Tinf; // Reference temperature [K]
551 
552  // Physical dimensions
553  amrex::Real tpi_Hpb; // PB height [m]
554  amrex::Real tpi_Lpb; // PB length [m]
555  amrex::Real tpi_Wpb; // PB width [m]
556  amrex::Real tpi_lref; // PB reference length [m]
557 
558  amrex::Real tpi_net_buoyant; // Perturbation net-zero calculation storage
559  amrex::Real tpi_pert_adjust; // Perturbation adjust for net-zero per cell adjustment
560 
561  // Perturbation data arrays
562  amrex::Vector<amrex::Real> pb_interval; // PB update time [s]
563  amrex::Vector<amrex::Real> pb_local_etime; // PB local elapsed time [s]
564  amrex::Vector<amrex::Real> pb_amp; // PB perturbation amplitude Ri:[K]
565  amrex::Vector<amrex::Real> pb_netZero; // PB array used for net zero sum calculation
566 
567  // Random number generation between range (used for interval calculation)
568  amrex::Real RandomReal (const amrex::Real min, const amrex::Real max)
569  {
570  amrex::Real r = (amrex::Real) rand() / (amrex::Real) RAND_MAX;
571  return min + r * (max - min);
572  }
573 
574 };
575 #endif
constexpr amrex::Real CONST_GRAV
Definition: ERF_Constants.H:21
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real pp(amrex::Real y)
Definition: ERF_MicrophysicsUtils.H:219
AMREX_FORCE_INLINE amrex::IntVect TileNoZ()
Definition: ERF_TileNoZ.H:11
Definition: ERF_TurbPertStruct.H:18
int tpi_layers
Definition: ERF_TurbPertStruct.H:541
void pseudoRandomPert(const int &boxIdx, const int &lev, const amrex::IndexType &m_ixtype)
Definition: ERF_TurbPertStruct.H:295
amrex::Real tpi_Tinf
Definition: ERF_TurbPertStruct.H:550
amrex::Real tpi_Wpb
Definition: ERF_TurbPertStruct.H:555
amrex::Real tpi_Lpb
Definition: ERF_TurbPertStruct.H:554
amrex::Real * get_pb_mag()
Definition: ERF_TurbPertStruct.H:504
amrex::Vector< amrex::Real > pb_local_etime
Definition: ERF_TurbPertStruct.H:563
amrex::Vector< amrex::BoxArray > pb_ba
Definition: ERF_TurbPertStruct.H:531
std::string pp_prefix
Definition: ERF_TurbPertStruct.H:526
amrex::Vector< amrex::Real > pb_amp
Definition: ERF_TurbPertStruct.H:564
void netZeroBuoyantAdjust(const int &boxIdx, const int &lev)
Definition: ERF_TurbPertStruct.H:357
void calc_tpi_update(const int lev, const amrex::Real dt, amrex::MultiFab &mf_xvel, amrex::MultiFab &mf_yvel, amrex::MultiFab &mf_cons)
Definition: ERF_TurbPertStruct.H:184
void calc_tpi_amp(const int &boxIdx, const amrex::Real &interval)
Definition: ERF_TurbPertStruct.H:269
amrex::MultiFab pb_cell
Definition: ERF_TurbPertStruct.H:536
amrex::Vector< amrex::Real > pb_netZero
Definition: ERF_TurbPertStruct.H:565
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)
Definition: ERF_TurbPertStruct.H:29
int tpi_offset
Definition: ERF_TurbPertStruct.H:542
amrex::Real tpi_Ti
Definition: ERF_TurbPertStruct.H:549
amrex::Real tpi_pert_adjust
Definition: ERF_TurbPertStruct.H:559
amrex::Vector< amrex::Real > pb_interval
Definition: ERF_TurbPertStruct.H:562
void debug(amrex::Real)
Definition: ERF_TurbPertStruct.H:508
amrex::Real * get_pb_netZero()
Definition: ERF_TurbPertStruct.H:505
amrex::Vector< amrex::Real > pb_mag
Definition: ERF_TurbPertStruct.H:532
amrex::Real tpi_lref
Definition: ERF_TurbPertStruct.H:556
amrex::Vector< int > tpi_boxDim
Definition: ERF_TurbPertStruct.H:544
amrex::Real tpi_nonDim
Definition: ERF_TurbPertStruct.H:548
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
int pt_type
Definition: ERF_TurbPertStruct.H:528
~TurbulentPerturbation()
Definition: ERF_TurbPertStruct.H:22
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
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)
Definition: ERF_TurbPertStruct.H:245
amrex::Real tpi_Hpb
Definition: ERF_TurbPertStruct.H:553
amrex::Vector< int > tpi_direction
Definition: ERF_TurbPertStruct.H:545