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.

33 {}

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.
415  {
416  for (int boxIdx = 0; boxIdx < pb_ba[lev].size(); boxIdx++) {
417  amrex::Box pbx = amrex::convert(pb_ba[lev][boxIdx], m_ixtype);
418  if (pt_type[lev] == 3) { pbx.setBig(2, pbx.bigEnd(2) - 1); } // prevent double counting after converting to k-faces
419  amrex::Box ubx = pbx & vbx;
420  if (ubx.ok()) {
421  if (comp == -1) { // vertical velocity perturbations
422  ParallelFor(ubx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
423  src_arr(i,j,k) += pert_cell(i,j,k);
424  });
425  } else {
426  ParallelFor(ubx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
427  src_arr(i,j,k,comp) += pert_cell(i,j,k);
428 
429  // For box region debug only
430  #ifdef INDEX_PERTURB
431  src_arr(i,j,k,comp) = (amrex::Real) (boxIdx + 5.);
432  #endif
433  });
434  }
435  }
436  }
437  }
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:757
amrex::Vector< int > pt_type
Integer perturbation type for each AMR level.
Definition: ERF_TurbPertStruct.H:754

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

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

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

◆ debug()

void TurbulentPerturbation::debug ( double  )
inline

Write perturbation debug information when debug output is enabled.

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

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

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

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

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.
491  {
492  for (amrex::MFIter mfi(pb_cell[lev],TileNoZ()); mfi.isValid(); ++mfi) {
493  amrex::Box vbx = mfi.validbox();
494  amrex::Box pbx = amrex::convert(pb_ba[lev][boxIdx], m_ixtype);
495  if (pt_type[lev] == 3) { pbx.setBig(2, pbx.bigEnd(2) - 1); } // prevent double counting after converting to k-faces
496  amrex::Box ubx = pbx & vbx;
497  if (ubx.ok()) {
498  amrex::Real amp_copy = pb_amp[lev][boxIdx];
499  const amrex::Array4<amrex::Real>& pert_cell = pb_cell[lev].array(mfi);
500 
501  if (pt_type[lev] == 2 || pt_type[lev] == 3) { // CPM
502  amrex::Real rand_number_const = RandomReal(-one, one);
503  ParallelFor(ubx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
504  pert_cell(i,j,k) = rand_number_const * amp_copy;
505  });
506  } else {
507  ParallelForRNG(ubx, [=] AMREX_GPU_DEVICE(int i, int j, int k, const amrex::RandomEngine& engine) noexcept {
508  amrex::Real rand_double = amrex::Random(engine);
509  pert_cell(i,j,k) = (rand_double*two - one) * amp_copy;
510  });
511  }
512  }
513  }
514  }
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) { const Real z_sfc=(use_eb) ? zero :fourth *(z_nd(i, j, klo)+z_nd(i+1, j, klo)+z_nd(i, j+1, klo)+z_nd(i+1, j+1, klo));const Real z_agl=z_cc(i, j, k) - z_sfc;state_pert(i, j, k, RhoKE_comp) *=amrex::max(std::pow(1 - amrex::min(z_agl/KE_decay_height, amrex::Real(1)), KE_decay_order), amrex::Real(1e-12));} } })
constexpr amrex::Real two
Definition: ERF_NumericalConstants.H:31

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.
807  {
808  amrex::Real r = (amrex::Real) rand() / (amrex::Real) RAND_MAX;
809  return min + r * (max - min);
810  }

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.
525  {
526 
527  for (amrex::MFIter mfi(pb_cell[lev],TileNoZ()); mfi.isValid(); ++mfi) {
528  amrex::Box vbx = mfi.validbox();
529  amrex::Box pbx = amrex::convert(pb_ba[lev][boxIdx], m_ixtype);
530  if (pt_type[lev] == 3) { pbx.setBig(2, pbx.bigEnd(2) - 1); } // prevent double counting after converting to k-faces
531  amrex::Box ubx = pbx & vbx;
532  if (ubx.ok()) {
533  const amrex::Array4<amrex::Real>& pert_cell = pb_cell[lev].array(mfi);
534  ParallelFor(ubx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
535  pert_cell(i,j,k) = zero;
536  });
537  }
538  }
539  }

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: