ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
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 
18 AMREX_ENUM(PerturbationType,
19  Source, Direct, CPM, None
20 );
21 
23 
24  public:
25 
27 
28  void init_tpi_type (const PerturbationType& pert_type)
29  {
30  if (pert_type == PerturbationType::Source) {
31  pt_type = 0;
32  } else if (pert_type == PerturbationType::Direct) {
33  pt_type = 1;
34  } else if (pert_type == PerturbationType::CPM) {
35  pt_type = 2;
36  }
37  AMREX_ALWAYS_ASSERT(pt_type >= 0);
38  }
39 
40  // Initializing Perturbation Region
41  // Currently only support perturbation in the x and y direction
42  // This portion of the code is only called once for initialization,
43  // therefore efficiency is not required. Dont spend too much time
44  // here other than to correctly setup the perturbation box region
45  void init_tpi (const int lev,
46  const amrex::IntVect& valid_box_lo, // start points of the valid box
47  const amrex::IntVect& valid_box_hi, // end points of the valid box
48  const amrex::GpuArray<amrex::Real,3> dx,
49  const amrex::BoxArray& ba,
50  const amrex::DistributionMapping& dm,
51  const int ngrow_state,
52  std::string pp_prefix,
53  const amrex::Vector<amrex::IntVect> refRatio,
54  const int max_level)
55 
56  {
57  // Initialization for some 0 dependent terms
58  tpi_Ti = 0.;
59  tpi_Tinf = 300.;
60 
61  ref_ratio = refRatio;
62 
63  amrex::ParmParse pp(pp_prefix);
64 
65  // Reading inputs, and placing assertion for the perturbation inflow to work
66  pp.queryarr("perturbation_box_dims",tpi_boxDim);
67  pp.queryarr("perturbation_direction",tpi_direction);
68  pp.query("perturbation_layers",tpi_layers);
69  pp.query("perturbation_offset",tpi_offset);
70 
71  pp.query("perturbation_nondimensional",tpi_nonDim);
72  pp.query("perturbation_T_infinity",tpi_Tinf);
73  pp.query("perturbation_T_intensity",tpi_Ti);
74  pp.query("perturbation_Ug",input_Ug);
75 
76  // Check variables message
77  if (tpi_layers < 0) { amrex::Abort("Please provide a valid perturbation layer value (ie. 3-5)"); }
78  if (tpi_offset < 0) { amrex::Abort("Please provide a valid inflow cell offset value for perturbation region (ie. 0-5)"); }
79  for (int i = 0; i < tpi_boxDim.size(); i++) {
80  if (tpi_boxDim[i] == 0) { amrex::Abort("Please provide valid dimensions for perturbation boxes."); }
81  }
82  if (tpi_nonDim < 0.) { amrex::Abort("Please provide a valid nondimensional number (ie. Ri = 0.042)"); }
83  if (input_Ug < 0.) { amrex::Abort("Please provide a valid geostrophic wind speed (ie. Ug = 10.0 m/s)"); }
84  if (tpi_Tinf < 0.) { amrex::Abort("Please provide a valid ambient temperature value (ie. T_0 = T_infty)"); }
85  if (tpi_Ti < 0.) { amrex::Abort("Please provide a valid temperature intensity value (ie. 0-1.0)"); }
86 
87  // Creating perturbation regions and initializing with generic size. Temporary size for now
88  amrex::Box lo_x_bx(amrex::IntVect(0), amrex::IntVect(1), amrex::IntVect(0));
89  amrex::Box hi_x_bx(amrex::IntVect(0), amrex::IntVect(1), amrex::IntVect(0));
90  amrex::Box lo_y_bx(amrex::IntVect(0), amrex::IntVect(1), amrex::IntVect(0));
91  amrex::Box hi_y_bx(amrex::IntVect(0), amrex::IntVect(1), amrex::IntVect(0));
92 
93  // Create a temporary box list to accumulate all the perturbation regions after box modification
94  amrex::BoxList tmp_bl;
95 
96  // boxSize for individual boxes
97  amrex::IntVect boxSize(tpi_boxDim[0],tpi_boxDim[1],tpi_boxDim[2]);
98 
99  // Starting logic to set the size of the perturbation region(s)
100  //amrex::PrintToFile("BoxPerturbationOutput") << "Setting perturbation region in:";
101  // ***** X-direction perturbation *****
102  if (tpi_direction[0]) { // West
103  lo_x_bx.setSmall(amrex::IntVect(valid_box_lo[0]+tpi_offset, valid_box_lo[1]+tpi_direction[1]*tpi_offset, valid_box_lo[2]));
104  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), valid_box_hi[2]));
105  //amrex::PrintToFile("BoxPerturbationOutput") << " West face";
106  }
107 
108  if (tpi_direction[3]) { // East
109  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, valid_box_lo[2]));
110  hi_x_bx.setBig(amrex::IntVect(valid_box_hi[0]-tpi_offset, valid_box_hi[1]-(tpi_direction[4]*tpi_offset), valid_box_hi[2]));
111  //amrex::PrintToFile("BoxPerturbationOutput") << " East face";
112  }
113 
114  // ***** Y-direction Perturbation *****
115  if (tpi_direction[1]) { // North
116  lo_y_bx.setSmall(amrex::IntVect(valid_box_lo[0]+tpi_direction[0]*tpi_offset, valid_box_lo[1]+tpi_offset, valid_box_lo[2]));
117  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, valid_box_hi[2]));
118  //amrex::PrintToFile("BoxPerturbationOutput") << " North face";
119  }
120 
121  if (tpi_direction[4]) { // South
122  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), valid_box_lo[2]));
123  hi_y_bx.setBig(amrex::IntVect(valid_box_hi[0]-tpi_direction[3]*tpi_offset, valid_box_hi[1]-tpi_offset, valid_box_hi[2]));
124  //amrex::PrintToFile("BoxPerturbationOutput") << " South face";
125  }
126 
127  if (tpi_direction[2] || tpi_direction[5]) { amrex::Abort("Currently not supporting z-direction flow perturbation"); }
128 
129  // Performing box union for intersecting perturbation regions to avoid overlapping sections (double counting at corners)
130  if (tpi_direction[0] && tpi_direction[1]) { // Reshaping South smallEnd
131  amrex::Box lo_x_lo_y_u = lo_x_bx & lo_y_bx;
132  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)));
133  }
134 
135  if (tpi_direction[3] && tpi_direction[1]) { // Reshaping South bigEnd
136  amrex::Box hi_x_lo_y_u = hi_x_bx & lo_y_bx;
137  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)));
138  }
139 
140  if (tpi_direction[0] && tpi_direction[4]) { // Reshaping North smallEnd
141  amrex::Box lo_x_hi_y_u = lo_x_bx & hi_y_bx;
142  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)));
143  }
144 
145  if (tpi_direction[3] && tpi_direction[4]) { // Reshaping North bigEnd
146  amrex::Box hi_x_hi_y_u = hi_x_bx & hi_y_bx;
147  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)));
148  }
149 
150  // Creating structure box array for conserved quantity
151  if (tpi_direction[0]) { tmp_bl.push_back(lo_x_bx); }
152  if (tpi_direction[1]) { tmp_bl.push_back(lo_y_bx); }
153  if (tpi_direction[3]) { tmp_bl.push_back(hi_x_bx); }
154  if (tpi_direction[4]) { tmp_bl.push_back(hi_y_bx); }
155  //amrex::PrintToFile("BoxPerturbationOutput") << "\nBoxList: " << tmp_bl << "\n";
156  amrex::BoxArray tmp_ba(tmp_bl);
157  tmp_ba.maxSize(boxSize);
158  pb_ba.push_back(tmp_ba);
159 
160  if (lev == 0) {
161  pb_mag.resize(max_level+1);
162  pb_dir.resize(max_level+1);
163  pb_netZero.resize(max_level+1);
164  pb_interval.resize(max_level+1);
165  pb_local_etime.resize(max_level+1);
166  pb_amp.resize(max_level+1);
167  pb_cell.resize(max_level+1);
168  tpi_Lpb.resize(max_level+1);
169  tpi_Wpb.resize(max_level+1);
170  tpi_Hpb.resize(max_level+1);
171  tpi_lref.resize(max_level+1);
172  }
173 
174  // Initializing mean magnitude and direction vectors
175  pb_mag[lev].resize(pb_ba[lev].size(), 0.);
176  pb_dir[lev].resize(pb_ba[lev].size(), 0.);
177  pb_netZero[lev].resize(pb_ba[lev].size(), 0.);
178 
179  // Set size of vector and initialize
180  pb_interval[lev].resize(pb_ba[lev].size(), -1.0);
181  pb_local_etime[lev].resize(pb_ba[lev].size(), 0.0);
182  pb_amp[lev].resize(pb_ba[lev].size(), 0.0);
183 
184  // Creating data array for perturbation amplitude storage
185  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...
186  pb_cell[lev].setVal(0.);
187 
188  // Computing perturbation reference length
189  tpi_Lpb[lev] = tpi_boxDim[0]*dx[0];
190  tpi_Wpb[lev] = tpi_boxDim[1]*dx[1];
191  tpi_Hpb[lev] = tpi_boxDim[2]*dx[2];
192  tpi_lref[lev] = sqrt(tpi_Lpb[lev]*tpi_Lpb[lev] + tpi_Wpb[lev]*tpi_Wpb[lev]);
193 
194  tpi_pert_adjust = 0.;
195  tpi_net_buoyant = 0.;
196 
197  /*
198  // Function check point message
199  amrex::PrintToFile("BoxPerturbationOutput") << "perturbation_box_dims: "
200  << tpi_boxDim[0] << " "
201  << tpi_boxDim[1] << " "
202  << tpi_boxDim[2] << "\n";
203  amrex::PrintToFile("BoxPerturbationOutput") << "perturbation_direction: "
204  << tpi_direction[0] << " "
205  << tpi_direction[1] << " "
206  << tpi_direction[2] << "\n\n";
207  amrex::PrintToFile("BoxPerturbationOutput") << "perturbation_layers: " << tpi_layers << "\n";
208  amrex::PrintToFile("BoxPerturbationOutput") << "perturbation_offset: " << tpi_offset << "\n\n";
209  amrex::PrintToFile("BoxPerturbationOutput") << "perturbation_nondimensional: " << tpi_nonDim << "\n";
210  amrex::PrintToFile("BoxPerturbationOutput") << "perturbation_T_infinity: " << tpi_Tinf << "\n";
211  amrex::PrintToFile("BoxPerturbationOutput") << "perturbation_T_intensity: " << tpi_Ti << "\n";
212  amrex::PrintToFile("BoxPerturbationOutput") << "Reference length per box = " << tpi_lref[lev] << "\n\n";
213  amrex::PrintToFile("BoxPerturbationOutput") << "Turbulent perturbation BoxArray:\n" << pb_ba[lev] << "\n";
214  */
215  }
216 
217 
218  // Perturbation update frequency check.
219  // This function will trigger calc_tpi_meanMag() and calc_tpi_amp(), when
220  // the local elapsed time is greater than the perturbation frequency.
221  // At each timestep a total buoyant force sanity check is made to ensure
222  // that the sum of the buoyant force introduced into the system is net-zero
223  void calc_tpi_update (const int lev,
224  const amrex::Real dt,
225  amrex::MultiFab& mf_xvel,
226  amrex::MultiFab& mf_yvel,
227  amrex::MultiFab& mf_cons)
228  {
229  // Resettubg the net buoyant force value
230  tpi_net_buoyant = 0.;
231 
232  // Setting random number generator for update interval
233  srand( (unsigned) time(NULL) );
234 
235  auto m_ixtype = mf_cons.boxArray().ixType(); // safety step
236 
237  // Seed the random generator at 1024UL for regression testing
238  int fix_random_seed = 0;
239  amrex::ParmParse pp("erf");
240  pp.query("fix_random_seed", fix_random_seed);
241  if (fix_random_seed) {
242  // We need this one for the ParalleForRNG used in calc_tpi
243  amrex::InitRandom(1024UL);
244 
245  // We need this one for the RandomReal below
246  srand(1024UL);
247  }
248 
249  for (int boxIdx = 0; boxIdx < pb_ba[lev].size(); boxIdx++) {
250 
251  bool update_box = true; // initialize flag
252  // Check if the local elapsed time is greater than the update interval and don't go into boxes the rank doesn't own
253  if (pt_type == 2) {
254  // for CPM, perturbation boxes/cells refresh after boxes have advected box width/length * num_layers (advective time scale)
255  update_box = ((pb_local_etime[lev][boxIdx] >= pb_interval[lev][boxIdx]*tpi_layers) && pb_local_etime[lev][boxIdx] != -1.0);
256  } else {
257  update_box = ( pb_local_etime[lev][boxIdx] >= pb_interval[lev][boxIdx] && pb_local_etime[lev][boxIdx] != -1.0 );
258  }
259  if ( update_box ) {
260 
261  // Compute mean velocity of each perturbation box
262  calc_tpi_meanMag_perBox(boxIdx, lev, mf_cons, mf_xvel, mf_yvel);
263 
264  // Only the rank owning the box will be able to access the storage location
265  // Done for parallelism to avoid Inf being stored in array
266  if (pb_mag[lev][boxIdx] !=0.) {
267  amrex::Real interval = 0.0;
268  if (pt_type == 2) {
269  // Wind direction correction for angled wind
270  amrex::Real wind_direction = pb_dir[lev][boxIdx];
271  if (wind_direction > PI / 4) { wind_direction = PI / 2 - wind_direction; }
272  // CPM only cares about the side length of a box, min call maintains flexibility.
273  interval = 1.0 / std::cos(wind_direction) * std::min(tpi_Lpb[lev], tpi_Wpb[lev]) / pb_mag[lev][boxIdx];
274  } else {
275  interval = tpi_lref[lev] / pb_mag[lev][boxIdx];
276  }
277  pb_interval[lev][boxIdx] = RandomReal(0.9*interval,1.1*interval); // 10% variation
278 
279  // Reset local elapsed time
280  pb_local_etime[lev][boxIdx] = 0.;
281  } else {
282  // this box is not on this rank, we shouldn't enter it again.
283  // Technically, all boxes are looped through on the very first step of a simulation and this is when this is set
284  pb_local_etime[lev][boxIdx] = -1.0;
285  }
286 
287  // Trigger amplitude calculation per perturbation box
288  calc_tpi_amp(lev, boxIdx, pb_interval[lev][boxIdx]);
289 
290  // Trigger random amplitude storage per cell within perturbation box
291  pseudoRandomPert(boxIdx, lev, m_ixtype);
292 
293  } else {
294  // set perturbation amplitudes to 0 for CPM. A little inefficient but leverages as much as existing code as possible.
295  if (pt_type == 2) { zero_amp(boxIdx, lev, m_ixtype); }
296 
297  // Increase by timestep of level 0 (but only if the box is owned by this rank)
298  if (pb_local_etime[lev][boxIdx] != -1.0) { pb_local_etime[lev][boxIdx] += dt; }
299  } // if
300 
301  if (pt_type < 2) { // box perturbation method only
302  // Per iteration operation of net-zero buoyant force check
303  if (pb_mag[lev][boxIdx] !=0.) netZeroBuoyantAdd(boxIdx, lev);
304  tpi_net_buoyant += pb_netZero[lev][boxIdx];
305  }
306  } // for
307 
308  if (pt_type < 2) { // box perturbation method only
309  // Normalizing the adjustment based on how many boxes there are
310  // the values within the array is already normalized by the number
311  // of cells within each box
313 
314  // Per iteration operation of net-zero buoyant force adjustment
315  for (int boxIdx = 0; boxIdx < pb_ba[lev].size(); boxIdx++) {
316  if (pb_mag[lev][boxIdx] !=0.) netZeroBuoyantAdjust(boxIdx, lev);
317  }
318  }
319  }
320 
321  // Applying perturbation amplitude onto source term (Umphrey and Senocak 2016)
322  // This function does per cell application within the box union. Random perturbation
323  // is assigned in calc_tpi_update.
324  void apply_tpi (const int& lev,
325  const amrex::Box& vbx, // box union from upper level
326  const int& comp, // Component to modify
327  const amrex::IndexType& m_ixtype, // IntVect type of src_arr
328  const amrex::Array4<amrex::Real>& src_arr, // Array to apply perturbation
329  const amrex::Array4<amrex::Real const>& pert_cell)
330  {
331  for (int boxIdx = 0; boxIdx < pb_ba[lev].size(); boxIdx++) {
332  amrex::Box pbx = amrex::convert(pb_ba[lev][boxIdx], m_ixtype);
333  amrex::Box ubx = pbx & vbx;
334  if (ubx.ok()) {
335  ParallelFor(ubx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
336  src_arr(i,j,k,comp) += pert_cell(i,j,k);
337 
338  // For box region debug only
339  #ifdef INDEX_PERTURB
340  src_arr(i,j,k,comp) = (amrex::Real) (boxIdx + 5.);
341  #endif
342  });
343  }
344  }
345  }
346 
347  // Perturbation amplitude calculation
348  void calc_tpi_amp (const int& lev,
349  const int& boxIdx,
350  const amrex::Real& interval)
351  {
352  pb_amp[lev][boxIdx] = 0.; // Safety step
353  if (pt_type == 2) { // CPM
354  amrex::Real cp = 1004; // specific heat of air [J/(kg K)]
355  amrex::Real Ec = 0.2; // Eckert number
356  pb_amp[lev][boxIdx] = (input_Ug * input_Ug) / (Ec * cp);
357  } else { // box perturbation
358  amrex::Real Um = pb_mag[lev][boxIdx];
359  amrex::Real beta = 1./tpi_Tinf; // Thermal expansion coefficient
360 
361  // Pseudo Random temperature the ignores scale when mechanically tripping turbulence
363  // get total refinement ratio on each level
364  int total_ref_ratio = 1;
365  for (int level = lev; level >= 1; level--) {
366  total_ref_ratio *= ref_ratio[level-1][2];
367  }
368  // calculation needs to be scale-aware since the formulation relies on the physical size of the box
369  if (tpi_Ti > 0.) g = (tpi_nonDim * Um * Um) / (tpi_Ti * tpi_Hpb[lev]) * 1 / total_ref_ratio;
370 
371  // Ma and Senocak (2023) Eq. 8, solving for delta phi
372  pb_amp[lev][boxIdx] = (tpi_nonDim * Um * Um) / (g * beta * tpi_Hpb[lev]) * 1 / total_ref_ratio;
373 
374  if (pt_type == 0) {
375  // Performing this step converts the perturbation proportionality into
376  // the forcing term
377  // Ma & Senocak (2023) Eq. 7
378  pb_amp[lev][boxIdx] /= interval;
379  }
380  }
381  }
382 
383  // Assigns pseudo-random (ie. white noise) perturbation to a storage cell, this
384  // value is then held constant for the duration of the update interval and assigned onto
385  // the source term
386  void pseudoRandomPert (const int& boxIdx,
387  const int& lev,
388  const amrex::IndexType& m_ixtype)
389  {
390  for (amrex::MFIter mfi(pb_cell[lev],TileNoZ()); mfi.isValid(); ++mfi) {
391  amrex::Box vbx = mfi.validbox();
392  amrex::Box pbx = amrex::convert(pb_ba[lev][boxIdx], m_ixtype);
393  amrex::Box ubx = pbx & vbx;
394  if (ubx.ok()) {
395  amrex::Real amp_copy = pb_amp[lev][boxIdx];
396  const amrex::Array4<amrex::Real>& pert_cell = pb_cell[lev].array(mfi);
397 
398  if (pt_type == 2) { // CPM
399  amrex::Real rand_number_const = RandomReal(-1.0, 1.0);
400  ParallelFor(ubx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
401  pert_cell(i,j,k) = rand_number_const * amp_copy;
402  });
403  } else {
404  ParallelForRNG(ubx, [=] AMREX_GPU_DEVICE(int i, int j, int k, const amrex::RandomEngine& engine) noexcept {
405  amrex::Real rand_double = amrex::Random(engine);
406  pert_cell(i,j,k) = (rand_double*2.0 - 1.0) * amp_copy;
407  });
408  }
409 
410  }
411  }
412  }
413 
414  // convert pert_cell value to be 0 so that perturbation amplitudes do not compound for CPM
415  void zero_amp (const int& boxIdx,
416  const int& lev,
417  const amrex::IndexType& m_ixtype)
418  {
419 
420  for (amrex::MFIter mfi(pb_cell[lev],TileNoZ()); mfi.isValid(); ++mfi) {
421  amrex::Box vbx = mfi.validbox();
422  amrex::Box pbx = amrex::convert(pb_ba[lev][boxIdx], m_ixtype);
423  amrex::Box ubx = pbx & vbx;
424  if (ubx.ok()) {
425  const amrex::Array4<amrex::Real>& pert_cell = pb_cell[lev].array(mfi);
426  ParallelFor(ubx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
427  pert_cell(i,j,k) = 0.0;
428  });
429  }
430  }
431  }
432 
433  // Checks for net-zero buoyant force introduction into the system
434  void netZeroBuoyantAdd (const int& boxIdx,
435  const int& lev)
436  {
437  // Creating local copy of PB box array and magnitude
438  const amrex::BoxArray m_pb_ba = pb_ba[lev];
439  amrex::Real* m_pb_netZero = pb_netZero[lev].data();
440 
441  // Create device array for summation
442  amrex::Vector<amrex::Real> avg_h(1,0.);
443  amrex::Gpu::DeviceVector<amrex::Real> avg_d(1,0.);
444  amrex::Real* avg = avg_d.data();
445 
446  // Iterates through the cells of each box and sum the white noise perturbation
447  for (amrex::MFIter mfi(pb_cell[lev], TileNoZ()) ; mfi.isValid(); ++mfi) {
448  const amrex::Box& vbx = mfi.validbox();
449  amrex::Box pbx = amrex::convert(m_pb_ba[boxIdx], vbx.ixType());
450  amrex::Box ubx = pbx & vbx;
451  if (ubx.ok()) {
452  const amrex::Array4<const amrex::Real>& pert_cell = pb_cell[lev].const_array(mfi);
453  amrex::Real norm = 1.0 / static_cast<amrex::Real>(ubx.numPts());
454  ParallelFor(amrex::Gpu::KernelInfo().setReduction(true), ubx, [=]
455  AMREX_GPU_DEVICE(int i, int j, int k, amrex::Gpu::Handler const& handler) noexcept {
456  amrex::Gpu::deviceReduceSum(&avg[0], pert_cell(i,j,k)*norm, handler);
457  });
458  amrex::Gpu::copy(amrex::Gpu::deviceToHost, avg_d.begin(), avg_d.end(), avg_h.begin());
459 
460  // Assigning onto storage array
461  m_pb_netZero[boxIdx] = avg_h[0];
462  }
463  }
464  }
465 
466  // If it's not a net-zero buoyant force, then adjust all cell by a normalized value
467  // to achieve this logic
468  void netZeroBuoyantAdjust (const int& boxIdx,
469  const int& lev)
470  {
471  // Creating local copy of PB box array and magnitude
472  const amrex::BoxArray m_pb_ba = pb_ba[lev];
473  for (amrex::MFIter mfi(pb_cell[lev], TileNoZ()) ; mfi.isValid(); ++mfi) {
474  const amrex::Box& vbx = mfi.validbox();
475  amrex::Box pbx = amrex::convert(m_pb_ba[boxIdx], vbx.ixType());
476  amrex::Box ubx = pbx & vbx;
477  if (ubx.ok()) {
478  const amrex::Real adjust = tpi_pert_adjust;
479  const amrex::Array4<amrex::Real>& pert_cell = pb_cell[lev].array(mfi);
480  ParallelFor(ubx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
481  pert_cell(i,j,k) -= adjust;
482  });
483  }
484  }
485  }
486 
487 // TODO: Test the difference between these two for Source term perturbation
488 #define USE_VOLUME_AVERAGE
489  // Perturbation box mean velocity magnitude calculation
490  // This is pulled into the structure to also utilize during runtime
491  void calc_tpi_meanMag_perBox (const int& boxIdx,
492  const int& lev,
493  amrex::MultiFab& mf_cons,
494  amrex::MultiFab& mf_xvel,
495  amrex::MultiFab& mf_yvel)
496 
497  {
498  // Creating local copy of PB box array and magnitude
499  const amrex::BoxArray m_pb_ba = pb_ba[lev];
500  amrex::Real* m_pb_mag = pb_mag[lev].data();
501  amrex::Real* m_pb_dir = pb_dir[lev].data();
502  m_pb_mag[boxIdx] = 0.; // Safety step
503  m_pb_dir[boxIdx] = 0.; // Safety step
504 
505  // Storage of averages per PB
506  // Index: 0=u (vol/slab_lo), 1=v (vol/slab_lo)
507  // 2=u (slab_hi), 3=v (slab_hi)
508  int n_avg = 4;
509  amrex::Vector<amrex::Real> avg_h(n_avg,0.);
510  amrex::Gpu::DeviceVector<amrex::Real> avg_d(n_avg,0.);
511  amrex::Real* avg = avg_d.data();
512 
513  // Averaging u & v components in single MFIter
514  for (amrex::MFIter mfi(mf_cons, TileNoZ()); mfi.isValid(); ++mfi) {
515 
516  // CC valid box (inherited from mf_cons)
517  const amrex::Box& vbx = mfi.validbox();
518 
519  // Box logic for u velocity
520  auto ixtype_u = mf_xvel.boxArray().ixType();
521  amrex::Box vbx_u = amrex::convert(vbx,ixtype_u);
522  amrex::Box pbx_u = amrex::convert(m_pb_ba[boxIdx], ixtype_u);
523  amrex::Box ubx_u = pbx_u & vbx_u;
524 
525  // Box logic for v velocity
526  auto ixtype_v = mf_yvel.boxArray().ixType();
527  amrex::Box vbx_v = amrex::convert(vbx,ixtype_v);
528  amrex::Box pbx_v = amrex::convert(m_pb_ba[boxIdx], ixtype_v);
529  amrex::Box ubx_v = pbx_v & vbx_v;
530 
531  // Operation over box union (U)
532  if (ubx_u.ok()) {
533  const amrex::Array4<const amrex::Real>& xvel_arry = mf_xvel.const_array(mfi);
534 
535  #ifdef USE_VOLUME_AVERAGE
536  amrex::Real norm = 1.0 / static_cast<amrex::Real>(ubx_u.numPts());
537  ParallelFor(amrex::Gpu::KernelInfo().setReduction(true), ubx_u, [=]
538  AMREX_GPU_DEVICE(int i, int j, int k, amrex::Gpu::Handler const& handler) noexcept {
539  amrex::Gpu::deviceReduceSum(&avg[0], xvel_arry(i,j,k)*norm, handler);
540  });
541  #endif // USE_VOLUME_AVERAGE
542 
543  #ifdef USE_SLAB_AVERAGE
544  amrex::Box ubxSlab_lo = makeSlab(ubx_u,2,ubx_u.smallEnd(2));
545  amrex::Box ubxSlab_hi = makeSlab(ubx_u,2,ubx_u.bigEnd(2));
546  amrex::Real norm_lo = 1.0 / static_cast<amrex::Real>(ubxSlab_lo.numPts());
547  amrex::Real norm_hi = 1.0 / static_cast<amrex::Real>(ubxSlab_hi.numPts());
548 
549  // Average u in the low slab
550  ParallelFor(amrex::Gpu::KernelInfo().setReduction(true), ubxSlab_lo, [=]
551  AMREX_GPU_DEVICE(int i, int j, int k, amrex::Gpu::Handler const& handler) noexcept {
552  amrex::Gpu::deviceReduceSum(&avg[0], xvel_arry(i,j,k)*norm_lo, handler);
553  });
554 
555  // Average u in the high slab
556  ParallelFor(amrex::Gpu::KernelInfo().setReduction(true), ubxSlab_hi, [=]
557  AMREX_GPU_DEVICE(int i, int j, int k, amrex::Gpu::Handler const& handler) noexcept {
558  amrex::Gpu::deviceReduceSum(&avg[2], xvel_arry(i,j,k)*norm_hi, handler);
559  });
560  #endif // USE_SLAB_AVERAGE
561  } // if
562 
563  // Operation over box union (V)
564  if (ubx_v.ok()) {
565  const amrex::Array4<const amrex::Real>& yvel_arry = mf_yvel.const_array(mfi);
566 
567  #ifdef USE_VOLUME_AVERAGE
568  amrex::Real norm = 1.0 / static_cast<amrex::Real>(ubx_v.numPts());
569  ParallelFor(amrex::Gpu::KernelInfo().setReduction(true), ubx_v, [=]
570  AMREX_GPU_DEVICE(int i, int j, int k, amrex::Gpu::Handler const& handler) noexcept {
571  amrex::Gpu::deviceReduceSum(&avg[1], yvel_arry(i,j,k)*norm, handler);
572  });
573  #endif // USE_VOLUME_AVERAGE
574 
575  #ifdef USE_SLAB_AVERAGE
576  amrex::Box ubxSlab_lo = makeSlab(ubx_v,2,ubx_v.smallEnd(2));
577  amrex::Box ubxSlab_hi = makeSlab(ubx_v,2,ubx_v.bigEnd(2));
578  amrex::Real norm_lo = 1.0 / static_cast<amrex::Real>(ubxSlab_lo.numPts());
579  amrex::Real norm_hi = 1.0 / static_cast<amrex::Real>(ubxSlab_hi.numPts());
580 
581  // Average v in the low slab
582  ParallelFor(amrex::Gpu::KernelInfo().setReduction(true), ubxSlab_lo, [=]
583  AMREX_GPU_DEVICE(int i, int j, int k, amrex::Gpu::Handler const& handler) noexcept {
584  amrex::Gpu::deviceReduceSum(&avg[1], yvel_arry(i,j,k)*norm_lo, handler);
585  });
586 
587  // Average v in the high slab
588  ParallelFor(amrex::Gpu::KernelInfo().setReduction(true), ubxSlab_hi, [=]
589  AMREX_GPU_DEVICE(int i, int j, int k, amrex::Gpu::Handler const& handler) noexcept {
590  amrex::Gpu::deviceReduceSum(&avg[3], yvel_arry(i,j,k)*norm_hi, handler);
591  });
592  #endif // USE_SLAB_AVERAGE
593  } // if
594  } // MFIter
595 
596  // Copy from device back to host
597  amrex::Gpu::copy(amrex::Gpu::deviceToHost, avg_d.begin(), avg_d.end(), avg_h.begin());
598 
599  // Computing the average magnitude within PB
600  #ifdef USE_VOLUME_AVERAGE
601  m_pb_mag[boxIdx] = sqrt(avg_h[0]*avg_h[0] + avg_h[1]*avg_h[1]);
602  m_pb_dir[boxIdx] = std::atan(std::abs(avg_h[0]) / std::abs(avg_h[1]+std::numeric_limits<amrex::Real>::epsilon()));
603  #endif
604 
605  #ifdef USE_SLAB_AVERAGE
606  m_pb_mag[boxIdx] = 0.5*( sqrt(avg_h[0]*avg_h[0] + avg_h[1]*avg_h[1])
607  + sqrt(avg_h[2]*avg_h[2] + avg_h[3]*avg_h[3]));
608  m_pb_dir[boxIdx] = std::atan(std::abs(0.5*(avg_h[0]+avg_h[2])) / std::abs(0.5*(avg_h[1]+avg_h[3])+std::numeric_limits<amrex::Real>::epsilon()));
609  #endif
610  }
611 
612  // Output debug message into a file
613  void debug (amrex::Real /*time*/)
614  {
615  /*
616  amrex::PrintToFile("BoxPerturbationOutput") << "#################### PB output at time = "
617  << time << " ####################\n";
618  amrex::PrintToFile("BoxPerturbationOutput") << " Using type: " << pt_type << "\n";
619  amrex::PrintToFile("BoxPerturbationOutput") << " Net: " << tpi_net_buoyant << " Adjust : " << tpi_pert_adjust << "\n";
620  for (int i = 0; i < pb_mag.size(); i++) {
621  amrex::PrintToFile("BoxPerturbationOutput") << "[" << i
622  << "] pb_Umag=" << pb_mag[i]
623  << " | pb_interval=" << pb_interval[i]
624  << " (" << pb_local_etime[i]
625  << ") | pb_amp=" << pb_amp[i] << "\n";
626  }
627  amrex::PrintToFile("BoxPerturbationOutput") << "\n";
628  */
629  }
630 
631  int pt_type = -1;
632 
633  // Public data members
634  amrex::Vector<amrex::BoxArray> pb_ba; // PB box array
635  amrex::Vector<amrex::Vector<amrex::Real>> pb_mag; // BP mean magnitude [m/s]
636  amrex::Vector<amrex::Vector<amrex::Real>> pb_dir; // BP mean direction [deg]
637 
638  // Perturbation amplitude cell storage
639  // This is after random assignment of equation (10) in Ma and Senocak 2023
640  amrex::Vector<amrex::MultiFab> pb_cell;
641 
642  private:
643 
644  // Private data members
645  int tpi_layers; // Number of layers of perturbation boxes
646  int tpi_offset; // Cells to offset the start of the perturbation region
647 
648  amrex::Vector<int> tpi_boxDim; // Dimensions of each perturbation box
649  amrex::Vector<int> tpi_direction; // Direction of the peturbation
650 
651  // Richardson Formulation
652  amrex::Real tpi_nonDim; // Richardson number
653  amrex::Real tpi_Ti; // Temperature intensity value
654  amrex::Real tpi_Tinf; // Reference temperature [K]
655 
656  // Physical dimensions
657  amrex::Vector<amrex::Real> tpi_Hpb; // PB height [m]
658  amrex::Vector<amrex::Real> tpi_Lpb; // PB length [m]
659  amrex::Vector<amrex::Real> tpi_Wpb; // PB width [m]
660  amrex::Vector<amrex::Real> tpi_lref; // PB reference length [m]
661 
662  amrex::Real tpi_net_buoyant; // Perturbation net-zero calculation storage
663  amrex::Real tpi_pert_adjust; // Perturbation adjust for net-zero per cell adjustment
664 
665  amrex::Vector<amrex::IntVect> ref_ratio; // ref_ratio for multilevel run
666  amrex::Real input_Ug; // input geostrophic wind speed to scale perturbations when using CPM
667 
668  // Perturbation data arrays
669  amrex::Vector<amrex::Vector<amrex::Real>> pb_interval; // PB update time [s]
670  amrex::Vector<amrex::Vector<amrex::Real>> pb_local_etime; // PB local elapsed time [s]
671  amrex::Vector<amrex::Vector<amrex::Real>> pb_amp; // PB perturbation amplitude Ri:[K]
672  amrex::Vector<amrex::Vector<amrex::Real>> pb_netZero; // PB array used for net zero sum calculation
673 
674  // Random number generation between range (used for interval calculation)
676  {
677  amrex::Real r = (amrex::Real) rand() / (amrex::Real) RAND_MAX;
678  return min + r * (max - min);
679  }
680 
681 };
682 #endif
constexpr amrex::Real PI
Definition: ERF_Constants.H:6
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:230
amrex::Real Real
Definition: ERF_ShocInterface.H:16
AMREX_FORCE_INLINE amrex::IntVect TileNoZ()
Definition: ERF_TileNoZ.H:11
AMREX_ENUM(PerturbationType, Source, Direct, CPM, None)
real(c_double), parameter g
Definition: ERF_module_model_constants.F90:19
real(c_double), parameter epsilon
Definition: ERF_module_model_constants.F90:12
real(c_double), parameter cp
Definition: ERF_module_model_constants.F90:22
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real avg_h(amrex::Real XXXm, amrex::Real XXXp)
Definition: ERF_EBUtils.H:10
Definition: ERF_TurbPertStruct.H:22
amrex::Vector< amrex::MultiFab > pb_cell
Definition: ERF_TurbPertStruct.H:640
int tpi_layers
Definition: ERF_TurbPertStruct.H:645
void pseudoRandomPert(const int &boxIdx, const int &lev, const amrex::IndexType &m_ixtype)
Definition: ERF_TurbPertStruct.H:386
amrex::Real tpi_Tinf
Definition: ERF_TurbPertStruct.H:654
amrex::Vector< amrex::Vector< amrex::Real > > pb_local_etime
Definition: ERF_TurbPertStruct.H:670
amrex::Vector< amrex::IntVect > ref_ratio
Definition: ERF_TurbPertStruct.H:665
void calc_tpi_amp(const int &lev, const int &boxIdx, const amrex::Real &interval)
Definition: ERF_TurbPertStruct.H:348
amrex::Real input_Ug
Definition: ERF_TurbPertStruct.H:666
amrex::Vector< amrex::BoxArray > pb_ba
Definition: ERF_TurbPertStruct.H:634
void netZeroBuoyantAdjust(const int &boxIdx, const int &lev)
Definition: ERF_TurbPertStruct.H:468
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:223
amrex::Vector< amrex::Vector< amrex::Real > > pb_mag
Definition: ERF_TurbPertStruct.H:635
amrex::Vector< amrex::Real > tpi_lref
Definition: ERF_TurbPertStruct.H:660
int tpi_offset
Definition: ERF_TurbPertStruct.H:646
amrex::Vector< amrex::Real > tpi_Hpb
Definition: ERF_TurbPertStruct.H:657
void init_tpi_type(const PerturbationType &pert_type)
Definition: ERF_TurbPertStruct.H:28
amrex::Vector< amrex::Vector< amrex::Real > > pb_netZero
Definition: ERF_TurbPertStruct.H:672
amrex::Real tpi_Ti
Definition: ERF_TurbPertStruct.H:653
amrex::Real tpi_pert_adjust
Definition: ERF_TurbPertStruct.H:663
void debug(amrex::Real)
Definition: ERF_TurbPertStruct.H:613
amrex::Vector< amrex::Vector< amrex::Real > > pb_interval
Definition: ERF_TurbPertStruct.H:669
amrex::Vector< amrex::Vector< amrex::Real > > pb_dir
Definition: ERF_TurbPertStruct.H:636
void zero_amp(const int &boxIdx, const int &lev, const amrex::IndexType &m_ixtype)
Definition: ERF_TurbPertStruct.H:415
amrex::Vector< int > tpi_boxDim
Definition: ERF_TurbPertStruct.H:648
amrex::Vector< amrex::Real > tpi_Wpb
Definition: ERF_TurbPertStruct.H:659
amrex::Real tpi_nonDim
Definition: ERF_TurbPertStruct.H:652
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:491
void init_tpi(const int lev, const amrex::IntVect &valid_box_lo, const amrex::IntVect &valid_box_hi, 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)
Definition: ERF_TurbPertStruct.H:45
void netZeroBuoyantAdd(const int &boxIdx, const int &lev)
Definition: ERF_TurbPertStruct.H:434
int pt_type
Definition: ERF_TurbPertStruct.H:631
~TurbulentPerturbation()
Definition: ERF_TurbPertStruct.H:26
amrex::Real RandomReal(const amrex::Real min, const amrex::Real max)
Definition: ERF_TurbPertStruct.H:675
amrex::Vector< amrex::Real > tpi_Lpb
Definition: ERF_TurbPertStruct.H:658
amrex::Real tpi_net_buoyant
Definition: ERF_TurbPertStruct.H:662
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:324
amrex::Vector< int > tpi_direction
Definition: ERF_TurbPertStruct.H:649
amrex::Vector< amrex::Vector< amrex::Real > > pb_amp
Definition: ERF_TurbPertStruct.H:671