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

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

#include <ERF_TurbPertStruct.H>

Collaboration diagram for TurbulentPerturbation:

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)
 Initialize turbulent perturbation regions and per-box storage. More...
 
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)
 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 netZeroBuoyantAdd (const int &boxIdx, const int &lev)
 Accumulate the mean perturbation contribution for net-zero buoyancy enforcement. More...
 
void netZeroBuoyantAdjust (const int &boxIdx, const int &lev)
 Adjust perturbation cells so the net buoyant forcing is zero. More...
 
void calc_tpi_meanMag_perBox (const int &boxIdx, const int &lev, amrex::MultiFab &mf_cons, amrex::MultiFab &mf_xvel, amrex::MultiFab &mf_yvel)
 Compute mean horizontal velocity magnitude and direction for one perturbation box. More...
 
void debug (double)
 Write perturbation debug information when debug output is enabled. 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 Member Functions

amrex::Real RandomReal (const amrex::Real min, const amrex::Real max)
 Return a pseudo-random real value in a closed interval. More...
 

Private Attributes

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::Realtpi_Hpb
 Perturbation-box height for each AMR level [m]. More...
 
amrex::Vector< amrex::Realtpi_Lpb
 Perturbation-box length for each AMR level [m]. More...
 
amrex::Vector< amrex::Realtpi_Wpb
 Perturbation-box width for each AMR level [m]. More...
 
amrex::Vector< amrex::Realtpi_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...
 

Detailed Description

Runtime state and operations for turbulent perturbation forcing.

Constructor & Destructor Documentation

◆ ~TurbulentPerturbation()

TurbulentPerturbation::~TurbulentPerturbation ( )
inline

Destroy turbulent perturbation state.

32 {}

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

Apply stored turbulent perturbations to a source or state array.

Parameters
levAMR level index.
vbxValid box over which perturbations may be applied.
compComponent to modify, or -1 for vertical velocity perturbations.
m_ixtypeIndex type of the destination array.
src_arrDestination array receiving perturbation increments.
pert_cellStored perturbation amplitudes.
414  {
415  for (int boxIdx = 0; boxIdx < pb_ba[lev].size(); boxIdx++) {
416  amrex::Box pbx = amrex::convert(pb_ba[lev][boxIdx], m_ixtype);
417  if (pt_type[lev] == 3) { pbx.setBig(2, pbx.bigEnd(2) - 1); } // prevent double counting after converting to k-faces
418  amrex::Box ubx = pbx & vbx;
419  if (ubx.ok()) {
420  if (comp == -1) { // vertical velocity perturbations
421  ParallelFor(ubx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
422  src_arr(i,j,k) += pert_cell(i,j,k);
423  });
424  } else {
425  ParallelFor(ubx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
426  src_arr(i,j,k,comp) += pert_cell(i,j,k);
427 
428  // For box region debug only
429  #ifdef INDEX_PERTURB
430  src_arr(i,j,k,comp) = (amrex::Real) (boxIdx + 5.);
431  #endif
432  });
433  }
434  }
435  }
436  }
ParallelFor(fab_box, [=] AMREX_GPU_DEVICE(int i, int j, int k) { qrcuten_arr(i, j, k)=Real(0);qscuten_arr(i, j, k)=Real(0);qicuten_arr(i, j, k)=Real(0);})
amrex::Real Real
Definition: ERF_ShocInterface.H:19
amrex::Vector< amrex::BoxArray > pb_ba
Perturbation-box BoxArray for each AMR level.
Definition: ERF_TurbPertStruct.H:756
amrex::Vector< int > pt_type
Integer perturbation type for each AMR level.
Definition: ERF_TurbPertStruct.H:753

Referenced by make_sources().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ calc_tpi_amp()

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

Compute the perturbation amplitude for one perturbation box.

Parameters
levAMR level index.
boxIdxPerturbation-box index on the level.
intervalCurrent perturbation update interval [s].
447  {
448  pb_amp[lev][boxIdx] = zero; // Safety step
449  if (pt_type[lev] == 2) { // CPM
450  amrex::Real cp = 1004; // specific heat of air [J/(kg K)]
451  amrex::Real Ec = amrex::Real(0.2); // Eckert number
452  pb_amp[lev][boxIdx] = (input_Ug * input_Ug) / (Ec * cp);
453  } else if (pt_type[lev] == 3) {
454  pb_amp[lev][boxIdx] = input_w_amp;
455  } else { // box perturbation
456  amrex::Real Um = pb_mag[lev][boxIdx];
457  amrex::Real beta = one/tpi_Tinf; // Thermal expansion coefficient
458 
459  // Pseudo Random temperature the ignores scale when mechanically tripping turbulence
461  // get total refinement ratio on each level
462  int total_ref_ratio = 1;
463  for (int level = lev; level >= 1; level--) {
464  total_ref_ratio *= ref_ratio[level-1][2];
465  }
466  // calculation needs to be scale-aware since the formulation relies on the physical size of the box
467  if (tpi_Ti > zero) g = (tpi_nonDim * Um * Um) / (tpi_Ti * tpi_Hpb[lev]) * 1 / total_ref_ratio;
468 
469  // Ma and Senocak (2023) Eq. 8, solving for delta phi
470  pb_amp[lev][boxIdx] = (tpi_nonDim * Um * Um) / (g * beta * tpi_Hpb[lev]) * 1 / total_ref_ratio;
471 
472  if (pt_type[lev] == 0) {
473  // Performing this step converts the perturbation proportionality into
474  // the forcing term
475  // Ma & Senocak (2023) Eq. 7
476  pb_amp[lev][boxIdx] /= static_cast<amrex::Real>(interval);
477  }
478  }
479  }
constexpr amrex::Real one
Definition: ERF_Constants.H:9
constexpr amrex::Real zero
Definition: ERF_Constants.H:8
constexpr amrex::Real CONST_GRAV
Definition: ERF_Constants.H:64
amrex::Real beta
Definition: ERF_InitCustomPert_DataAssimilation_ISV.H:10
real(c_double), parameter g
Definition: ERF_module_model_constants.F90:19
real(c_double), parameter cp
Definition: ERF_module_model_constants.F90:22
amrex::Real tpi_Tinf
Reference temperature used by the perturbation formulation [K].
Definition: ERF_TurbPertStruct.H:776
amrex::Real input_w_amp
Input vertical velocity perturbation amplitude.
Definition: ERF_TurbPertStruct.H:789
amrex::Vector< amrex::Vector< amrex::Real > > pb_mag
Mean velocity magnitude for each perturbation box [m/s].
Definition: ERF_TurbPertStruct.H:757
amrex::Real input_Ug
Input geostrophic wind speed used to scale CPM perturbations.
Definition: ERF_TurbPertStruct.H:788
amrex::Vector< amrex::IntVect > ref_ratio
Refinement ratios used by multilevel perturbation scaling.
Definition: ERF_TurbPertStruct.H:787
amrex::Real tpi_Ti
Temperature intensity used by the perturbation formulation.
Definition: ERF_TurbPertStruct.H:775
amrex::Vector< amrex::Vector< amrex::Real > > pb_amp
Perturbation amplitude for each perturbation box.
Definition: ERF_TurbPertStruct.H:796
amrex::Real tpi_nonDim
Nondimensional number used by the perturbation formulation.
Definition: ERF_TurbPertStruct.H:774
amrex::Vector< amrex::Real > tpi_Hpb
Perturbation-box height for each AMR level [m].
Definition: ERF_TurbPertStruct.H:779

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

Compute mean horizontal velocity magnitude and direction for one perturbation box.

Parameters
boxIdxPerturbation-box index on the level.
levAMR level index.
mf_consConserved state used for iteration layout.
mf_xvelX-velocity field.
mf_yvelY-velocity field.
617  {
618  // Creating local copy of PB box array and magnitude
619  const amrex::BoxArray m_pb_ba = pb_ba[lev];
620  amrex::Real* m_pb_mag = pb_mag[lev].data();
621  amrex::Real* m_pb_dir = pb_dir[lev].data();
622  m_pb_mag[boxIdx] = zero; // Safety step
623  m_pb_dir[boxIdx] = zero; // Safety step
624 
625  // Storage of averages per PB
626  // Index: 0=u (vol/slab_lo), 1=v (vol/slab_lo)
627  // 2=u (slab_hi), 3=v (slab_hi)
628  int n_avg = 4;
629  amrex::Vector<amrex::Real> avg_h(n_avg,zero);
630  amrex::Gpu::DeviceVector<amrex::Real> avg_d(n_avg,zero);
631  amrex::Real* avg = avg_d.data();
632 
633  // Averaging u & v components in single MFIter
634  for (amrex::MFIter mfi(mf_cons, TileNoZ()); mfi.isValid(); ++mfi) {
635 
636  // CC valid box (inherited from mf_cons)
637  const amrex::Box& vbx = mfi.validbox();
638 
639  // Box logic for u velocity
640  auto ixtype_u = mf_xvel.boxArray().ixType();
641  amrex::Box vbx_u = amrex::convert(vbx,ixtype_u);
642  amrex::Box pbx_u = amrex::convert(m_pb_ba[boxIdx], ixtype_u);
643  amrex::Box ubx_u = pbx_u & vbx_u;
644 
645  // Box logic for v velocity
646  auto ixtype_v = mf_yvel.boxArray().ixType();
647  amrex::Box vbx_v = amrex::convert(vbx,ixtype_v);
648  amrex::Box pbx_v = amrex::convert(m_pb_ba[boxIdx], ixtype_v);
649  amrex::Box ubx_v = pbx_v & vbx_v;
650 
651  // Operation over box union (U)
652  if (ubx_u.ok()) {
653  const amrex::Array4<const amrex::Real>& xvel_arry = mf_xvel.const_array(mfi);
654 
655  #ifdef USE_VOLUME_AVERAGE
656  amrex::Real norm = one / static_cast<amrex::Real>(ubx_u.numPts());
657  ParallelFor(amrex::Gpu::KernelInfo().setReduction(true), ubx_u, [=]
658  AMREX_GPU_DEVICE(int i, int j, int k, amrex::Gpu::Handler const& handler) noexcept {
659  amrex::Gpu::deviceReduceSum(&avg[0], xvel_arry(i,j,k)*norm, handler);
660  });
661  #endif // USE_VOLUME_AVERAGE
662 
663  #ifdef USE_SLAB_AVERAGE
664  amrex::Box ubxSlab_lo = makeSlab(ubx_u,2,ubx_u.smallEnd(2));
665  amrex::Box ubxSlab_hi = makeSlab(ubx_u,2,ubx_u.bigEnd(2));
666  amrex::Real norm_lo = one / static_cast<amrex::Real>(ubxSlab_lo.numPts());
667  amrex::Real norm_hi = one / static_cast<amrex::Real>(ubxSlab_hi.numPts());
668 
669  // Average u in the low slab
670  ParallelFor(amrex::Gpu::KernelInfo().setReduction(true), ubxSlab_lo, [=]
671  AMREX_GPU_DEVICE(int i, int j, int k, amrex::Gpu::Handler const& handler) noexcept {
672  amrex::Gpu::deviceReduceSum(&avg[0], xvel_arry(i,j,k)*norm_lo, handler);
673  });
674 
675  // Average u in the high slab
676  ParallelFor(amrex::Gpu::KernelInfo().setReduction(true), ubxSlab_hi, [=]
677  AMREX_GPU_DEVICE(int i, int j, int k, amrex::Gpu::Handler const& handler) noexcept {
678  amrex::Gpu::deviceReduceSum(&avg[2], xvel_arry(i,j,k)*norm_hi, handler);
679  });
680  #endif // USE_SLAB_AVERAGE
681  } // if
682 
683  // Operation over box union (V)
684  if (ubx_v.ok()) {
685  const amrex::Array4<const amrex::Real>& yvel_arry = mf_yvel.const_array(mfi);
686 
687  #ifdef USE_VOLUME_AVERAGE
688  amrex::Real norm = one / static_cast<amrex::Real>(ubx_v.numPts());
689  ParallelFor(amrex::Gpu::KernelInfo().setReduction(true), ubx_v, [=]
690  AMREX_GPU_DEVICE(int i, int j, int k, amrex::Gpu::Handler const& handler) noexcept {
691  amrex::Gpu::deviceReduceSum(&avg[1], yvel_arry(i,j,k)*norm, handler);
692  });
693  #endif // USE_VOLUME_AVERAGE
694 
695  #ifdef USE_SLAB_AVERAGE
696  amrex::Box ubxSlab_lo = makeSlab(ubx_v,2,ubx_v.smallEnd(2));
697  amrex::Box ubxSlab_hi = makeSlab(ubx_v,2,ubx_v.bigEnd(2));
698  amrex::Real norm_lo = one / static_cast<amrex::Real>(ubxSlab_lo.numPts());
699  amrex::Real norm_hi = one / static_cast<amrex::Real>(ubxSlab_hi.numPts());
700 
701  // Average v in the low slab
702  ParallelFor(amrex::Gpu::KernelInfo().setReduction(true), ubxSlab_lo, [=]
703  AMREX_GPU_DEVICE(int i, int j, int k, amrex::Gpu::Handler const& handler) noexcept {
704  amrex::Gpu::deviceReduceSum(&avg[1], yvel_arry(i,j,k)*norm_lo, handler);
705  });
706 
707  // Average v in the high slab
708  ParallelFor(amrex::Gpu::KernelInfo().setReduction(true), ubxSlab_hi, [=]
709  AMREX_GPU_DEVICE(int i, int j, int k, amrex::Gpu::Handler const& handler) noexcept {
710  amrex::Gpu::deviceReduceSum(&avg[3], yvel_arry(i,j,k)*norm_hi, handler);
711  });
712  #endif // USE_SLAB_AVERAGE
713  } // if
714  } // MFIter
715 
716  // Copy from device back to host
717  amrex::Gpu::copy(amrex::Gpu::deviceToHost, avg_d.begin(), avg_d.end(), avg_h.begin());
718 
719  // Computing the average magnitude within PB
720  #ifdef USE_VOLUME_AVERAGE
721  m_pb_mag[boxIdx] = std::sqrt(avg_h[0]*avg_h[0] + avg_h[1]*avg_h[1]);
722  m_pb_dir[boxIdx] = std::atan(std::abs(avg_h[0]) / std::abs(avg_h[1]+std::numeric_limits<amrex::Real>::epsilon()));
723  #endif
724 
725  #ifdef USE_SLAB_AVERAGE
726  m_pb_mag[boxIdx] = myhalf*( std::sqrt(avg_h[0]*avg_h[0] + avg_h[1]*avg_h[1])
727  + std::sqrt(avg_h[2]*avg_h[2] + avg_h[3]*avg_h[3]));
728  m_pb_dir[boxIdx] = std::atan(std::abs(myhalf*(avg_h[0]+avg_h[2])) / std::abs(myhalf*(avg_h[1]+avg_h[3])+std::numeric_limits<amrex::Real>::epsilon()));
729  #endif
730  }
constexpr amrex::Real myhalf
Definition: ERF_Constants.H:13
AMREX_FORCE_INLINE amrex::IntVect TileNoZ()
Definition: ERF_TileNoZ.H:11
real(c_double), parameter epsilon
Definition: ERF_module_model_constants.F90:12
amrex::Vector< amrex::Vector< amrex::Real > > pb_dir
Mean velocity direction for each perturbation box.
Definition: ERF_TurbPertStruct.H:758

Referenced by calc_tpi_update().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ calc_tpi_update()

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

Update perturbation amplitudes and intervals when each box is ready.

Parameters
levAMR level index.
dtTime step for elapsed-time accumulation.
mf_xvelX-velocity field used to compute box mean velocities.
mf_yvelY-velocity field used to compute box mean velocities.
mf_consConserved state used for iteration layout.
304  {
305  // Resetting the net buoyant force value
307 
308  // Setting random number generator for update interval
309  srand( (unsigned) time(NULL) );
310 
311  auto m_ixtype = mf_cons.boxArray().ixType(); // safety step
312  if (pt_type[lev] == 3) { m_ixtype = amrex::IndexType(amrex::IntVect(0,0,1)); }
313 
314  // Seed the random generator at 1024UL for regression testing
315  int fix_random_seed = 0;
316  amrex::ParmParse pp("erf");
317  pp.queryAdd("fix_random_seed", fix_random_seed);
318  if (fix_random_seed) {
319  // We need this one for the ParalleForRNG used in calc_tpi
320  amrex::InitRandom(1024UL, amrex::ParallelDescriptor::NProcs(), 1024UL);
321 
322  // We need this one for the RandomReal below
323  srand(1024UL);
324  }
325 
326  for (int boxIdx = 0; boxIdx < pb_ba[lev].size(); boxIdx++) {
327 
328  bool update_box = true; // initialize flag
329  // Check if the local elapsed time is greater than the update interval and don't go into boxes the rank doesn't own
330  if (pt_type[lev] == 2 || pt_type[lev] == 3) {
331  // for CPM, perturbation boxes/cells refresh after boxes have advected box width/length * num_layers (advective time scale)
332  update_box = ((pb_local_etime[lev][boxIdx] >= pb_interval[lev][boxIdx]*tpi_layers) && pb_local_etime[lev][boxIdx] != -one);
333  } else {
334  update_box = ( pb_local_etime[lev][boxIdx] >= pb_interval[lev][boxIdx] && pb_local_etime[lev][boxIdx] != -one );
335  }
336  if ( update_box ) {
337 
338  // Compute mean velocity of each perturbation box
339  calc_tpi_meanMag_perBox(boxIdx, lev, mf_cons, mf_xvel, mf_yvel);
340 
341  // Only the rank owning the box will be able to access the storage location
342  // Done for parallelism to avoid Inf being stored in array
343  if (pb_mag[lev][boxIdx] !=zero) {
344  double interval = 0.0;
345  if (pt_type[lev] == 2 || pt_type[lev] == 3) {
346  // Wind direction correction for angled wind
347  amrex::Real wind_direction = pb_dir[lev][boxIdx];
348  if (wind_direction > PI / 4) { wind_direction = PI / 2 - wind_direction; }
349  // CPM only cares about the side length of a box, min call maintains flexibility.
350  interval = one / std::cos(wind_direction) * std::min(tpi_Lpb[lev], tpi_Wpb[lev]) / pb_mag[lev][boxIdx];
351  } else {
352  interval = tpi_lref[lev] / pb_mag[lev][boxIdx];
353  }
354  pb_interval[lev][boxIdx] = RandomReal(static_cast<amrex::Real>(0.9*interval),
355  static_cast<amrex::Real>(1.1*interval)); // 10% variation
356 
357  // Reset local elapsed time
358  pb_local_etime[lev][boxIdx] = zero;
359  } else {
360  // this box is not on this rank, we shouldn't enter it again.
361  // Technically, all boxes are looped through on the very first step of a simulation and this is when this is set
362  pb_local_etime[lev][boxIdx] = -one;
363  }
364 
365  // Trigger amplitude calculation per perturbation box
366  calc_tpi_amp(lev, boxIdx, pb_interval[lev][boxIdx]);
367 
368  // Trigger random amplitude storage per cell within perturbation box
369  pseudoRandomPert(boxIdx, lev, m_ixtype);
370 
371  } else {
372  // set perturbation amplitudes to 0 for CPM. A little inefficient but leverages as much as existing code as possible.
373  if (pt_type[lev] == 2 || pt_type[lev] == 3) { zero_amp(boxIdx, lev, m_ixtype); }
374 
375  // Increase by timestep of level 0 (but only if the box is owned by this rank)
376  if (pb_local_etime[lev][boxIdx] != -one) { pb_local_etime[lev][boxIdx] += dt; }
377  } // if
378 
379  if (pt_type[lev] < 2) { // box perturbation method only
380  // Per iteration operation of net-zero buoyant force check
381  if (pb_mag[lev][boxIdx] !=zero) netZeroBuoyantAdd(boxIdx, lev);
382  tpi_net_buoyant += pb_netZero[lev][boxIdx];
383  }
384  } // for
385 
386  if (pt_type[lev] < 2) { // box perturbation method only
387  // Normalizing the adjustment based on how many boxes there are
388  // the values within the array is already normalized by the number
389  // of cells within each box
391 
392  // Per iteration operation of net-zero buoyant force adjustment
393  for (int boxIdx = 0; boxIdx < pb_ba[lev].size(); boxIdx++) {
394  if (pb_mag[lev][boxIdx] !=zero) netZeroBuoyantAdjust(boxIdx, lev);
395  }
396  }
397  }
constexpr amrex::Real PI
Definition: ERF_Constants.H:42
ParmParse pp("prob")
int fix_random_seed
Definition: ERF_InitCustomPert_ABL.H:5
int tpi_layers
Number of layers of perturbation boxes.
Definition: ERF_TurbPertStruct.H:767
amrex::Vector< amrex::Real > tpi_Wpb
Perturbation-box width for each AMR level [m].
Definition: ERF_TurbPertStruct.H:781
void calc_tpi_amp(const int &lev, const int &boxIdx, const double &interval)
Compute the perturbation amplitude for one perturbation box.
Definition: ERF_TurbPertStruct.H:444
void pseudoRandomPert(const int &boxIdx, const int &lev, const amrex::IndexType &m_ixtype)
Assign pseudo-random perturbations to cells in one perturbation box.
Definition: ERF_TurbPertStruct.H:487
amrex::Real tpi_pert_adjust
Per-cell perturbation adjustment used for net-zero buoyancy.
Definition: ERF_TurbPertStruct.H:785
amrex::Vector< amrex::Real > tpi_Lpb
Perturbation-box length for each AMR level [m].
Definition: ERF_TurbPertStruct.H:780
amrex::Vector< amrex::Vector< amrex::Real > > pb_netZero
Per-box storage used for net-zero buoyancy calculation.
Definition: ERF_TurbPertStruct.H:797
void netZeroBuoyantAdjust(const int &boxIdx, const int &lev)
Adjust perturbation cells so the net buoyant forcing is zero.
Definition: ERF_TurbPertStruct.H:582
amrex::Real tpi_net_buoyant
Accumulated net buoyant perturbation used for correction.
Definition: ERF_TurbPertStruct.H:784
amrex::Vector< amrex::Vector< double > > pb_interval
Perturbation update interval for each box [s].
Definition: ERF_TurbPertStruct.H:794
amrex::Vector< amrex::Vector< double > > pb_local_etime
Local elapsed time for each perturbation box [s].
Definition: ERF_TurbPertStruct.H:795
amrex::Vector< amrex::Real > tpi_lref
Perturbation-box reference length for each AMR level [m].
Definition: ERF_TurbPertStruct.H:782
void zero_amp(const int &boxIdx, const int &lev, const amrex::IndexType &m_ixtype)
Reset CPM perturbation amplitudes in one perturbation box to zero.
Definition: ERF_TurbPertStruct.H:521
void calc_tpi_meanMag_perBox(const int &boxIdx, const int &lev, amrex::MultiFab &mf_cons, amrex::MultiFab &mf_xvel, amrex::MultiFab &mf_yvel)
Compute mean horizontal velocity magnitude and direction for one perturbation box.
Definition: ERF_TurbPertStruct.H:611
void netZeroBuoyantAdd(const int &boxIdx, const int &lev)
Accumulate the mean perturbation contribution for net-zero buoyancy enforcement.
Definition: ERF_TurbPertStruct.H:545
amrex::Real RandomReal(const amrex::Real min, const amrex::Real max)
Return a pseudo-random real value in a closed interval.
Definition: ERF_TurbPertStruct.H:805
Here is the call graph for this function:

◆ debug()

void TurbulentPerturbation::debug ( double  )
inline

Write perturbation debug information when debug output is enabled.

736  {
737  /*
738  amrex::PrintToFile("BoxPerturbationOutput") << "#################### PB output at time = "
739  << time << " ####################\n";
740  amrex::PrintToFile("BoxPerturbationOutput") << " Using type: " << pt_type << "\n";
741  amrex::PrintToFile("BoxPerturbationOutput") << " Net: " << tpi_net_buoyant << " Adjust : " << tpi_pert_adjust << "\n";
742  for (int i = 0; i < pb_mag.size(); i++) {
743  amrex::PrintToFile("BoxPerturbationOutput") << "[" << i
744  << "] pb_Umag=" << pb_mag[i]
745  << " | pb_interval=" << pb_interval[i]
746  << " (" << pb_local_etime[i]
747  << ") | pb_amp=" << pb_amp[i] << "\n";
748  }
749  amrex::PrintToFile("BoxPerturbationOutput") << "\n";
750  */
751  }

◆ init_tpi()

void TurbulentPerturbation::init_tpi ( const int  lev,
const amrex::Vector< amrex::BoxArray > &  subdomains_lev,
const amrex::GpuArray< amrex::Real, 3 >  dx,
const amrex::BoxArray &  ba,
const amrex::DistributionMapping &  dm,
const int  ngrow_state,
std::string  pp_prefix,
const amrex::Vector< amrex::IntVect >  refRatio,
const int  max_level 
)
inline

Initialize turbulent perturbation regions and per-box storage.

Parameters
levAMR level index.
subdomains_levRectangular level subdomains used to define perturbation regions.
dxCell spacing for this level.
baState BoxArray used for perturbation-cell storage.
dmDistributionMapping for perturbation-cell storage.
ngrow_stateNumber of ghost cells for perturbation-cell storage.
pp_prefixParmParse prefix for the ERF input namespace.
refRatioRefinement ratios between AMR levels.
max_levelMaximum AMR level configured for the run.
83  {
84  ref_ratio = refRatio;
85 
86  amrex::ParmParse pp(pp_prefix);
87 
88  // Reading inputs, and placing assertion for the perturbation inflow to work
89  pp.getarr("perturbation_box_dims",tpi_boxDim);
90  pp.getarr("perturbation_direction",tpi_direction);
91  pp.get("perturbation_layers",tpi_layers);
92  pp.get("perturbation_offset",tpi_offset);
93 
94  tpi_nonDim = zero;
95  pp.queryAdd("perturbation_nondimensional",tpi_nonDim);
96 
97  tpi_Tinf = amrex::Real(300.);
98  pp.queryAdd("perturbation_T_infinity",tpi_Tinf);
99 
100  tpi_Ti = zero;
101  pp.queryAdd("perturbation_T_intensity",tpi_Ti);
102 
103  input_Ug = zero;
104  pp.queryAdd("perturbation_Ug",input_Ug);
105 
106  input_w_amp = zero;
107  pp.queryAdd("perturbation_w_amp",input_w_amp);
108 
109  // Negative sentinels rather than the queryAdd return values: k = 0 is a legal
110  // bound, so the natural default cannot double as "unset", and the return value
111  // only reports whether the key existed before the call. init_tpi runs once per
112  // level, so testing the return would make every level above 0 see "specified"
113  // with the value level 0 planted, skipping the valid_box_lo/hi defaults below.
114  perturbation_klo = -1;
115  pp.queryAdd("perturbation_klo", perturbation_klo);
116  const bool have_klo = (perturbation_klo >= 0);
117 
118  perturbation_khi = -1;
119  pp.queryAdd("perturbation_khi", perturbation_khi);
120  const bool have_khi = (perturbation_khi >= 0);
121 
122  // Check variables message
123  if (tpi_offset < 0) { amrex::Abort("Please provide a valid inflow cell offset value for perturbation region (ie. 0-5)"); }
124  if (tpi_layers < 0) { amrex::Abort("Please provide a valid perturbation layer value (ie. 3-5)"); }
125  if (tpi_nonDim < zero) { amrex::Abort("Please provide a valid nondimensional number (ie. Ri = amrex::Real(0.042))"); }
126  for (int i = 0; i < tpi_boxDim.size(); i++) {
127  if (tpi_boxDim[i] < 3) { amrex::Abort("Please provide valid dimensions for perturbation boxes."); }
128  }
129  if (input_Ug < zero) { amrex::Abort("Please provide a valid geostrophic wind speed (ie. Ug = amrex::Real(10.0) m/s)"); }
130  if (tpi_Tinf < zero) { amrex::Abort("Please provide a valid ambient temperature value (ie. T_0 = T_infty)"); }
131  if (tpi_Ti < zero) { amrex::Abort("Please provide a valid temperature intensity value (ie. 0-one)"); }
132 
133  // Create a temporary box list to accumulate all the perturbation regions after box modification
134  amrex::BoxList tmp_bl;
135 
136  // boxSize for individual boxes
137  amrex::IntVect boxSize(tpi_boxDim[0],tpi_boxDim[1],tpi_boxDim[2]);
138 
139  if (tpi_direction[2] || tpi_direction[5]) { amrex::Abort("Currently not supporting z-direction flow perturbation"); }
140 
141  for (int isub = 0; isub < subdomains_lev.size(); ++isub) {
142  const amrex::BoxArray& subdomain = subdomains_lev[isub];
143  amrex::Box subdomain_box(subdomain.minimalBox());
144 
145  if (subdomain_box.numPts() != subdomain.numPts()) {
146  amrex::Abort("Turbulent perturbations require rectangular subdomains. "
147  "Level " + std::to_string(lev) +
148  ", subdomain " + std::to_string(isub) +
149  " is not a rectangular region fully covered by grids.");
150  }
151 
152  const amrex::IntVect& valid_box_lo = subdomain_box.smallEnd();
153  const amrex::IntVect& valid_box_hi = subdomain_box.bigEnd();
154 
155  // default perturbation region to be the entire z extent if the user does not specify bounds
156  if (!have_klo) { perturbation_klo = valid_box_lo[2]; }
157  if (!have_khi) { perturbation_khi = valid_box_hi[2]; }
158 
159  // Creating perturbation regions and initializing with generic size.
160  amrex::Box lo_x_bx(amrex::IntVect(0), amrex::IntVect(1), amrex::IntVect(0));
161  amrex::Box hi_x_bx(amrex::IntVect(0), amrex::IntVect(1), amrex::IntVect(0));
162  amrex::Box lo_y_bx(amrex::IntVect(0), amrex::IntVect(1), amrex::IntVect(0));
163  amrex::Box hi_y_bx(amrex::IntVect(0), amrex::IntVect(1), amrex::IntVect(0));
164 
165  // Starting logic to set the size of the perturbation region(s)
166  //amrex::PrintToFile("BoxPerturbationOutput") << "Setting perturbation region in:";
167  // ***** X-direction perturbation *****
168  if (tpi_direction[0]) { // West
169  lo_x_bx.setSmall(amrex::IntVect(valid_box_lo[0]+tpi_offset, valid_box_lo[1]+tpi_direction[1]*tpi_offset, perturbation_klo));
170  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));
171  amrex::PrintToFile("BoxPerturbationOutput") << " West face";
172  }
173 
174  if (tpi_direction[3]) { // East
175  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));
176  hi_x_bx.setBig(amrex::IntVect(valid_box_hi[0]-tpi_offset, valid_box_hi[1]-(tpi_direction[4]*tpi_offset), perturbation_khi));
177  amrex::PrintToFile("BoxPerturbationOutput") << " East face";
178  }
179 
180  // ***** Y-direction Perturbation *****
181  if (tpi_direction[1]) { // North
182  lo_y_bx.setSmall(amrex::IntVect(valid_box_lo[0]+tpi_direction[0]*tpi_offset, valid_box_lo[1]+tpi_offset, perturbation_klo));
183  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));
184  amrex::PrintToFile("BoxPerturbationOutput") << " North face";
185  }
186 
187  if (tpi_direction[4]) { // South
188  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));
189  hi_y_bx.setBig(amrex::IntVect(valid_box_hi[0]-tpi_direction[3]*tpi_offset, valid_box_hi[1]-tpi_offset, perturbation_khi));
190  amrex::PrintToFile("BoxPerturbationOutput") << " South face";
191  }
192 
193  // Performing box union for intersecting perturbation regions to avoid overlapping sections (double counting at corners)
194  if (tpi_direction[0] && tpi_direction[1]) { // Reshaping South smallEnd
195  amrex::Box lo_x_lo_y_u = lo_x_bx & lo_y_bx;
196  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)));
197  }
198 
199  if (tpi_direction[3] && tpi_direction[1]) { // Reshaping South bigEnd
200  amrex::Box hi_x_lo_y_u = hi_x_bx & lo_y_bx;
201  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)));
202  }
203 
204  if (tpi_direction[0] && tpi_direction[4]) { // Reshaping North smallEnd
205  amrex::Box lo_x_hi_y_u = lo_x_bx & hi_y_bx;
206  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)));
207  }
208 
209  if (tpi_direction[3] && tpi_direction[4]) { // Reshaping North bigEnd
210  amrex::Box hi_x_hi_y_u = hi_x_bx & hi_y_bx;
211  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)));
212  }
213 
214  // Creating structure box array for conserved quantity
215  if (tpi_direction[0]) { tmp_bl.push_back(lo_x_bx); }
216  if (tpi_direction[1]) { tmp_bl.push_back(lo_y_bx); }
217  if (tpi_direction[3]) { tmp_bl.push_back(hi_x_bx); }
218  if (tpi_direction[4]) { tmp_bl.push_back(hi_y_bx); }
219  }
220 
221  //amrex::PrintToFile("BoxPerturbationOutput") << "\nBoxList: " << tmp_bl << "\n";
222  amrex::BoxArray tmp_ba(tmp_bl);
223  tmp_ba.maxSize(boxSize);
224 
225  const int num_levels = max_level + 1;
226  if (pb_ba.size() < num_levels) {
227  pb_ba.resize(num_levels);
228  pb_mag.resize(num_levels);
229  pb_dir.resize(num_levels);
230  pb_netZero.resize(num_levels);
231  pb_interval.resize(num_levels);
232  pb_local_etime.resize(num_levels);
233  pb_amp.resize(num_levels);
234  pb_cell.resize(num_levels);
235  tpi_Lpb.resize(num_levels);
236  tpi_Wpb.resize(num_levels);
237  tpi_Hpb.resize(num_levels);
238  tpi_lref.resize(num_levels);
239  }
240 
241  pb_ba[lev] = tmp_ba;
242 
243  // Initializing mean magnitude and direction vectors
244  pb_mag[lev].resize(pb_ba[lev].size(), zero);
245  pb_dir[lev].resize(pb_ba[lev].size(), zero);
246  pb_netZero[lev].resize(pb_ba[lev].size(), zero);
247 
248  // Set size of vector and initialize
249  pb_interval[lev].resize(pb_ba[lev].size(), -one);
250  pb_local_etime[lev].resize(pb_ba[lev].size(), zero);
251  pb_amp[lev].resize(pb_ba[lev].size(), zero);
252 
253  // Creating data array for perturbation amplitude storage
254  if (pt_type[lev] == 3) { // CPM_W converts to the k-face ba.
255  pb_cell[lev].define(convert(ba, amrex::IntVect(0,0,1)), dm, 1, ngrow_state);
256  } else {
257  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...
258  }
259  pb_cell[lev].setVal(0.);
260 
261  // Computing perturbation reference length
262  tpi_Lpb[lev] = tpi_boxDim[0]*dx[0];
263  tpi_Wpb[lev] = tpi_boxDim[1]*dx[1];
264  tpi_Hpb[lev] = tpi_boxDim[2]*dx[2];
265  tpi_lref[lev] = std::sqrt(tpi_Lpb[lev]*tpi_Lpb[lev] + tpi_Wpb[lev]*tpi_Wpb[lev]);
266 
269 
270  /*
271  // Function check point message
272  amrex::PrintToFile("BoxPerturbationOutput") << "perturbation_box_dims: "
273  << tpi_boxDim[0] << " "
274  << tpi_boxDim[1] << " "
275  << tpi_boxDim[2] << "\n";
276  amrex::PrintToFile("BoxPerturbationOutput") << "perturbation_direction: "
277  << tpi_direction[0] << " "
278  << tpi_direction[1] << " "
279  << tpi_direction[2] << "\n\n";
280  amrex::PrintToFile("BoxPerturbationOutput") << "perturbation_layers: " << tpi_layers << "\n";
281  amrex::PrintToFile("BoxPerturbationOutput") << "perturbation_offset: " << tpi_offset << "\n\n";
282  amrex::PrintToFile("BoxPerturbationOutput") << "perturbation_nondimensional: " << tpi_nonDim << "\n";
283  amrex::PrintToFile("BoxPerturbationOutput") << "perturbation_T_infinity: " << tpi_Tinf << "\n";
284  amrex::PrintToFile("BoxPerturbationOutput") << "perturbation_T_intensity: " << tpi_Ti << "\n";
285  amrex::PrintToFile("BoxPerturbationOutput") << "Reference length per box = " << tpi_lref[lev] << "\n\n";
286  amrex::PrintToFile("BoxPerturbationOutput") << "Turbulent perturbation BoxArray:\n" << pb_ba[lev] << "\n";
287  */
288  }
const Real dx
Definition: ERF_InitCustomPert_ABL.H:44
integer, private isub
Definition: ERF_module_mp_morr_two_moment.F90:164
amrex::Vector< int > tpi_direction
Boundary directions where perturbations are applied.
Definition: ERF_TurbPertStruct.H:771
int perturbation_klo
Lower vertical index bound for perturbations.
Definition: ERF_TurbPertStruct.H:790
int tpi_offset
Cell offset for the start of the perturbation region.
Definition: ERF_TurbPertStruct.H:768
amrex::Vector< int > tpi_boxDim
Dimensions of each perturbation box.
Definition: ERF_TurbPertStruct.H:770
amrex::Vector< amrex::MultiFab > pb_cell
Per-cell perturbation amplitude storage.
Definition: ERF_TurbPertStruct.H:762
int perturbation_khi
Upper vertical index bound for perturbations.
Definition: ERF_TurbPertStruct.H:791
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
levAMR level index.
pert_typePerturbation type selected for this level.
max_levelMaximum AMR level configured for the run.
43  {
44  if (pt_type.size() < max_level + 1) {
45  pt_type.resize(max_level + 1, -1);
46  }
47 
48  if (pert_type == PerturbationType::Source) {
49  pt_type[lev] = 0;
50  } else if (pert_type == PerturbationType::Direct) {
51  pt_type[lev] = 1;
52  } else if (pert_type == PerturbationType::CPM) {
53  pt_type[lev] = 2;
54  } else if (pert_type == PerturbationType::CPM_W) {
55  pt_type[lev] = 3;
56  } else {
57  pt_type[lev] = -1;
58  }
59  }

◆ netZeroBuoyantAdd()

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

Accumulate the mean perturbation contribution for net-zero buoyancy enforcement.

Parameters
boxIdxPerturbation-box index on the level.
levAMR level index.
547  {
548  // Creating local copy of PB box array and magnitude
549  const amrex::BoxArray m_pb_ba = pb_ba[lev];
550  amrex::Real* m_pb_netZero = pb_netZero[lev].data();
551 
552  // Create device array for summation
553  amrex::Vector<amrex::Real> avg_h(1,zero);
554  amrex::Gpu::DeviceVector<amrex::Real> avg_d(1,zero);
555  amrex::Real* avg = avg_d.data();
556 
557  // Iterates through the cells of each box and sum the white noise perturbation
558  for (amrex::MFIter mfi(pb_cell[lev], TileNoZ()) ; mfi.isValid(); ++mfi) {
559  const amrex::Box& vbx = mfi.validbox();
560  amrex::Box pbx = amrex::convert(m_pb_ba[boxIdx], vbx.ixType());
561  amrex::Box ubx = pbx & vbx;
562  if (ubx.ok()) {
563  const amrex::Array4<const amrex::Real>& pert_cell = pb_cell[lev].const_array(mfi);
564  amrex::Real norm = one / static_cast<amrex::Real>(ubx.numPts());
565  ParallelFor(amrex::Gpu::KernelInfo().setReduction(true), ubx, [=]
566  AMREX_GPU_DEVICE(int i, int j, int k, amrex::Gpu::Handler const& handler) noexcept {
567  amrex::Gpu::deviceReduceSum(&avg[0], pert_cell(i,j,k)*norm, handler);
568  });
569  amrex::Gpu::copy(amrex::Gpu::deviceToHost, avg_d.begin(), avg_d.end(), avg_h.begin());
570 
571  // Assigning onto storage array
572  m_pb_netZero[boxIdx] = avg_h[0];
573  }
574  }
575  }

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

Adjust perturbation cells so the net buoyant forcing is zero.

Parameters
boxIdxPerturbation-box index on the level.
levAMR level index.
584  {
585  // Creating local copy of PB box array and magnitude
586  const amrex::BoxArray m_pb_ba = pb_ba[lev];
587  for (amrex::MFIter mfi(pb_cell[lev], TileNoZ()) ; mfi.isValid(); ++mfi) {
588  const amrex::Box& vbx = mfi.validbox();
589  amrex::Box pbx = amrex::convert(m_pb_ba[boxIdx], vbx.ixType());
590  amrex::Box ubx = pbx & vbx;
591  if (ubx.ok()) {
592  const amrex::Real adjust = tpi_pert_adjust;
593  const amrex::Array4<amrex::Real>& pert_cell = pb_cell[lev].array(mfi);
594  ParallelFor(ubx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
595  pert_cell(i,j,k) -= adjust;
596  });
597  }
598  }
599  }

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
boxIdxPerturbation-box index on the level.
levAMR level index.
m_ixtypeIndex type used to align perturbation storage.
490  {
491  for (amrex::MFIter mfi(pb_cell[lev],TileNoZ()); mfi.isValid(); ++mfi) {
492  amrex::Box vbx = mfi.validbox();
493  amrex::Box pbx = amrex::convert(pb_ba[lev][boxIdx], m_ixtype);
494  if (pt_type[lev] == 3) { pbx.setBig(2, pbx.bigEnd(2) - 1); } // prevent double counting after converting to k-faces
495  amrex::Box ubx = pbx & vbx;
496  if (ubx.ok()) {
497  amrex::Real amp_copy = pb_amp[lev][boxIdx];
498  const amrex::Array4<amrex::Real>& pert_cell = pb_cell[lev].array(mfi);
499 
500  if (pt_type[lev] == 2 || pt_type[lev] == 3) { // CPM
501  amrex::Real rand_number_const = RandomReal(-one, one);
502  ParallelFor(ubx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
503  pert_cell(i,j,k) = rand_number_const * amp_copy;
504  });
505  } else {
506  ParallelForRNG(ubx, [=] AMREX_GPU_DEVICE(int i, int j, int k, const amrex::RandomEngine& engine) noexcept {
507  amrex::Real rand_double = amrex::Random(engine);
508  pert_cell(i,j,k) = (rand_double*two - one) * amp_copy;
509  });
510  }
511  }
512  }
513  }
constexpr amrex::Real two
Definition: ERF_Constants.H:10
ParallelForRNG(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k, const amrex::RandomEngine &engine) noexcept { const Real x=prob_lo_x+(i+myhalf) *dx;const Real y=prob_lo_y+(j+myhalf) *dy;const Real z=z_cc(i, j, k);const Real r=std::sqrt((x-xc) *(x-xc)+(y-yc) *(y-yc)+(z-zc) *(z-zc));if((z<=pert_ref_height) &&(T_0_Pert_Mag !=amrex::Real(0))) { Real rand_double;if(deterministic_ic_pert) { rand_double=erf_hash_rng::hash_uniform(i, j, k, temperature_hash_comp, lev, ic_pert_seed);} else { rand_double=amrex::Random(engine);} state_pert(i, j, k, RhoTheta_comp)=(rand_double *amrex::Real(2) - amrex::Real(1)) *T_0_Pert_Mag;if(!pert_rhotheta) { state_pert(i, j, k, RhoTheta_comp) *=r_hse(i, j, k);} } if(add_sinusoidal_T) { const Real zl=z/pert_ref_height;const Real damp=std::exp(-myhalf *zl *zl);const Real rho=state(i, j, k, Rho_comp);const Real rhotheta=state(i, j, k, RhoTheta_comp);const Real qv=(use_moisture) ? state(i, j, k, RhoQ1_comp)/rho :amrex::Real(0);const Real Told=getTgivenRandRTh(rho, rhotheta, qv);const Real P=getPgivenRTh(rhotheta, qv);const Real Tnew=Told+tfac *damp *z *std::cos(cval *(x - xc));const Real theta_new=getThgivenTandP(Tnew, P, rdOcp);const Real rho_new=getRhogivenThetaPress(theta_new, P, rdOcp, qv);state_pert(i, j, k, Rho_comp)+=rho_new - rho;if(use_moisture) { state_pert(i, j, k, RhoQ1_comp)+=(rho_new - rho) *qv;} } state_pert(i, j, k, RhoScalar_comp)=A_0 *std::exp(-amrex::Real(10.) *r *r);if(state_pert.nComp() > RhoKE_comp) { if(rhoKE_0 > 0) { state_pert(i, j, k, RhoKE_comp)=rhoKE_0;} else { state_pert(i, j, k, RhoKE_comp)=r_hse(i, j, k) *KE_0;} if(KE_decay_height > 0) { state_pert(i, j, k, RhoKE_comp) *=amrex::max(std::pow(1 - amrex::min(z/KE_decay_height, amrex::Real(1)), KE_decay_order), Real(1e-12));} } })

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

Return a pseudo-random real value in a closed interval.

Parameters
minLower bound of the interval.
maxUpper bound of the interval.
Returns
Pseudo-random value between min and max.
806  {
807  amrex::Real r = (amrex::Real) rand() / (amrex::Real) RAND_MAX;
808  return min + r * (max - min);
809  }

Referenced by calc_tpi_update(), and pseudoRandomPert().

Here is the caller graph for this function:

◆ zero_amp()

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

Reset CPM perturbation amplitudes in one perturbation box to zero.

Parameters
boxIdxPerturbation-box index on the level.
levAMR level index.
m_ixtypeIndex type used to align perturbation storage.
524  {
525 
526  for (amrex::MFIter mfi(pb_cell[lev],TileNoZ()); mfi.isValid(); ++mfi) {
527  amrex::Box vbx = mfi.validbox();
528  amrex::Box pbx = amrex::convert(pb_ba[lev][boxIdx], m_ixtype);
529  if (pt_type[lev] == 3) { pbx.setBig(2, pbx.bigEnd(2) - 1); } // prevent double counting after converting to k-faces
530  amrex::Box ubx = pbx & vbx;
531  if (ubx.ok()) {
532  const amrex::Array4<amrex::Real>& pert_cell = pb_cell[lev].array(mfi);
533  ParallelFor(ubx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
534  pert_cell(i,j,k) = zero;
535  });
536  }
537  }
538  }

Referenced by calc_tpi_update().

Here is the call graph for this function:
Here is the caller graph for this function:

Member Data Documentation

◆ input_Ug

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

◆ pb_cell

TurbulentPerturbation::pb_cell

Per-cell perturbation amplitude storage.

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

◆ pb_dir

TurbulentPerturbation::pb_dir

Mean velocity direction for each perturbation box.

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

◆ pb_interval

TurbulentPerturbation::pb_interval
private

Perturbation update interval for each box [s].

Referenced by calc_tpi_update(), and init_tpi().

◆ pb_local_etime

TurbulentPerturbation::pb_local_etime
private

Local elapsed time for each perturbation box [s].

Referenced by calc_tpi_update(), 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_perBox(), calc_tpi_update(), and init_tpi().

◆ pb_netZero

TurbulentPerturbation::pb_netZero
private

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

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

◆ 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(), init_tpi(), init_tpi_type(), pseudoRandomPert(), 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 init_tpi().

◆ tpi_direction

TurbulentPerturbation::tpi_direction
private

Boundary directions where perturbations are applied.

Referenced by init_tpi().

◆ 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 calc_tpi_update(), and init_tpi().

◆ tpi_Lpb

TurbulentPerturbation::tpi_Lpb
private

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

Referenced by calc_tpi_update(), and init_tpi().

◆ 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 calc_tpi_update(), and init_tpi().

◆ 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 init_tpi().

◆ tpi_pert_adjust

TurbulentPerturbation::tpi_pert_adjust
private

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

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

◆ 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(), and init_tpi().


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