ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
Functions | Variables
ERF_InteriorGhostCells.cpp File Reference
#include <ERF_Utils.H>
Include dependency graph for ERF_InteriorGhostCells.cpp:

Functions

void realbdy_interior_bxs_xy (const Box &bx, const Box &domain, const int &width, Box &bx_xlo, Box &bx_xhi, Box &bx_ylo, Box &bx_yhi, const IntVect &ng_vect, const bool get_int_ng)
 
void realbdy_compute_interior_ghost_rhs (const Real &time, const Real &delta_t, const Real &start_bdy_time, const Real &final_bdy_time, const Real &bdy_time_interval, const Real &nudge_factor, int width, bool do_upwind, const Geometry &geom, Vector< MultiFab > &S_rhs, Vector< MultiFab > &S_cur_data, Vector< Vector< FArrayBox >> &bdy_data_xlo, Vector< Vector< FArrayBox >> &bdy_data_xhi, Vector< Vector< FArrayBox >> &bdy_data_ylo, Vector< Vector< FArrayBox >> &bdy_data_yhi, std::unique_ptr< ReadBndryPlanes > &m_r2d)
 
void fine_compute_interior_ghost_rhs (const Real &time, const Real &delta_t, const int &width, const int &set_width, const Geometry &geom, ERFFillPatcher *FPr_c, ERFFillPatcher *FPr_u, ERFFillPatcher *FPr_v, ERFFillPatcher *FPr_w, Vector< BCRec > &domain_bcs_type, Vector< MultiFab > &S_rhs_f, Vector< MultiFab > &S_data_f)
 

Variables

PhysBCFunctNoOp void_bc
 

Function Documentation

◆ fine_compute_interior_ghost_rhs()

void fine_compute_interior_ghost_rhs ( const Real &  time,
const Real &  delta_t,
const int &  width,
const int &  set_width,
const Geometry &  geom,
ERFFillPatcher *  FPr_c,
ERFFillPatcher *  FPr_u,
ERFFillPatcher *  FPr_v,
ERFFillPatcher *  FPr_w,
Vector< BCRec > &  domain_bcs_type,
Vector< MultiFab > &  S_rhs_f,
Vector< MultiFab > &  S_data_f 
)

Compute the RHS in the fine relaxation zone

Parameters
[in]timecurrent (elapsed) time
[in]delta_ttimestep
[in]widthnumber of cells in (relaxation+specified) zone
[in]set_widthnumber of cells in (specified) zone
[in]FPr_ccons fine patch container
[in]FPr_uuvel fine patch container
[in]FPr_vvvel fine patch container
[in]FPr_wwvel fine patch container
[in]boxes_at_levelboxes at current level
[in]domain_bcs_typeboundary condition types
[out]S_rhsRHS to be computed here
[in]S_datacurrent value of the solution
517 {
518  BL_PROFILE_REGION("fine_compute_interior_ghost_RHS()");
519 
520  // Relaxation constants
521  Real F1 = 1./(10.*delta_t);
522  Real F2 = 1./(50.*delta_t);
523 
524  // Vector of MFs to hold data (dm differs w/ fine patch)
525  Vector<MultiFab> fmf_p_v;
526 
527  // Loop over the variables
528  for (int ivar_idx = 0; ivar_idx < IntVars::NumTypes; ++ivar_idx)
529  {
530  // Fine mfs
531  MultiFab& fmf = S_data_f[ivar_idx];
532  MultiFab& rhs = S_rhs_f [ivar_idx];
533 
534  // NOTE: These temporary MFs and copy operations are horrible
535  // for memory usage and efficiency. However, we need to
536  // have access to ghost cells in the cons array to convert
537  // from primitive u/v/w to momentum. Furthermore, the BA
538  // for the fine patches in ERFFillPatcher don't match the
539  // BA for the data/RHS. For this reason, the data is copied
540  // to a vector of MFs (with ghost cells) so the BAs match
541  // the BA of data/RHS and we have access to rho to convert
542  // prim to conserved.
543 
544  // Temp MF on box (distribution map differs w/ fine patch)
545  int num_var = fmf.nComp();
546  fmf_p_v.emplace_back(fmf.boxArray(), fmf.DistributionMap(), num_var, fmf.nGrowVect());
547  MultiFab& fmf_p = fmf_p_v[ivar_idx];
548  MultiFab::Copy(fmf_p,fmf, 0, 0, num_var, fmf.nGrowVect());
549 
550  // Integer mask MF
551  int set_mask_val;
552  int relax_mask_val;
553  iMultiFab* mask;
554 
555  // Fill fine patch on interior halo region
556  //==========================================================
557  if (ivar_idx == IntVars::cons)
558  {
559  FPr_c->FillRelax(fmf_p, time, void_bc, domain_bcs_type);
560  mask = FPr_c->GetMask();
561  set_mask_val = FPr_c->GetSetMaskVal();
562  relax_mask_val = FPr_c->GetRelaxMaskVal();
563  }
564  else if (ivar_idx == IntVars::xmom)
565  {
566  FPr_u->FillRelax(fmf_p, time, void_bc, domain_bcs_type);
567  mask = FPr_u->GetMask();
568  set_mask_val = FPr_u->GetSetMaskVal();
569  relax_mask_val = FPr_u->GetRelaxMaskVal();
570 
571 #ifdef _OPENMP
572 #pragma omp parallel if (Gpu::notInLaunchRegion())
573 #endif
574  for ( MFIter mfi(fmf_p,TilingIfNotGPU()); mfi.isValid(); ++mfi)
575  {
576  Box tbx = mfi.tilebox();
577 
578  const Array4<Real>& prim_arr = fmf_p.array(mfi);
579  const Array4<const Real>& rho_arr = fmf_p_v[0].const_array(mfi);
580  const Array4<const int>& mask_arr = mask->const_array(mfi);
581 
582  ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
583  {
584  if (mask_arr(i,j,k) == relax_mask_val) {
585  Real rho_interp = 0.5 * ( rho_arr(i-1,j,k) + rho_arr(i,j,k) );
586  prim_arr(i,j,k) *= rho_interp;
587  }
588  });
589  } // mfi
590  }
591  else if (ivar_idx == IntVars::ymom)
592  {
593  FPr_v->FillRelax(fmf_p, time, void_bc, domain_bcs_type);
594  mask = FPr_v->GetMask();
595  set_mask_val = FPr_v->GetSetMaskVal();
596  relax_mask_val = FPr_v->GetRelaxMaskVal();
597 
598 #ifdef _OPENMP
599 #pragma omp parallel if (Gpu::notInLaunchRegion())
600 #endif
601  for ( MFIter mfi(fmf_p,TilingIfNotGPU()); mfi.isValid(); ++mfi)
602  {
603  Box tbx = mfi.tilebox();
604 
605  const Array4<Real>& prim_arr = fmf_p.array(mfi);
606  const Array4<const Real>& rho_arr = fmf_p_v[0].const_array(mfi);
607  const Array4<const int>& mask_arr = mask->const_array(mfi);
608 
609  ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
610  {
611  if (mask_arr(i,j,k) == relax_mask_val) {
612  Real rho_interp = 0.5 * ( rho_arr(i,j-1,k) + rho_arr(i,j,k) );
613  prim_arr(i,j,k) *= rho_interp;
614  }
615  });
616  } // mfi
617  }
618  else if (ivar_idx == IntVars::zmom)
619  {
620  FPr_w->FillRelax(fmf_p, time, void_bc, domain_bcs_type);
621  mask = FPr_w->GetMask();
622  set_mask_val = FPr_w->GetSetMaskVal();
623  relax_mask_val = FPr_w->GetRelaxMaskVal();
624 
625 #ifdef _OPENMP
626 #pragma omp parallel if (Gpu::notInLaunchRegion())
627 #endif
628  for ( MFIter mfi(fmf_p,TilingIfNotGPU()); mfi.isValid(); ++mfi)
629  {
630  Box tbx = mfi.tilebox();
631 
632  const Array4<Real>& prim_arr = fmf_p.array(mfi);
633  const Array4<const Real>& rho_arr = fmf_p_v[0].const_array(mfi);
634  const Array4<const int>& mask_arr = mask->const_array(mfi);
635 
636  ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
637  {
638  if (mask_arr(i,j,k) == relax_mask_val) {
639  Real rho_interp = 0.5 * ( rho_arr(i,j,k-1) + rho_arr(i,j,k) );
640  prim_arr(i,j,k) *= rho_interp;
641  }
642  });
643  } // mfi
644  } else {
645  Abort("Dont recognize this variable type in fine_compute_interior_ghost_RHS");
646  }
647 
648 
649  // Zero RHS in set region
650  //==========================================================
651 #ifdef _OPENMP
652 #pragma omp parallel if (Gpu::notInLaunchRegion())
653 #endif
654  for ( MFIter mfi(rhs,TilingIfNotGPU()); mfi.isValid(); ++mfi)
655  {
656  Box tbx = mfi.tilebox();
657  const Array4<Real>& rhs_arr = rhs.array(mfi);
658  const Array4<const int>& mask_arr = mask->const_array(mfi);
659 
660  ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
661  {
662  if (mask_arr(i,j,k) == set_mask_val) {
663  rhs_arr(i,j,k) = 0.0;
664  }
665  });
666  } // mfi
667 
668  // For Laplacian stencil
669  rhs.FillBoundary(geom.periodicity());
670 
671 
672  // Compute RHS in relaxation region
673  //==========================================================
674 #ifdef _OPENMP
675 #pragma omp parallel if (Gpu::notInLaunchRegion())
676 #endif
677  for ( MFIter mfi(fmf_p,TilingIfNotGPU()); mfi.isValid(); ++mfi)
678  {
679  Box tbx = mfi.tilebox();
680  const Array4<Real>& rhs_arr = rhs.array(mfi);
681  const Array4<const Real>& fine_arr = fmf_p.const_array(mfi);
682  const Array4<const Real>& data_arr = fmf.const_array(mfi);
683  const Array4<const int>& mask_arr = mask->const_array(mfi);
684 
685  Box vbx = mfi.validbox();
686  const auto& vbx_lo = lbound(vbx);
687  const auto& vbx_hi = ubound(vbx);
688 
689  int icomp = 0;
690 
691  int Spec_z = set_width;
692  int Relax_z = width - Spec_z;
693  Real num = Real(Spec_z + Relax_z);
694  Real denom = Real(Relax_z - 1);
695  ParallelFor(tbx, num_var, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
696  {
697  if (mask_arr(i,j,k) == relax_mask_val) {
698 
699  // Indices
700  Real n_ind(-1); // Set to -1 to quiet compiler warning
701  int ii(width-1); int jj(width-1);
702  bool near_x_lo_wall(false); bool near_x_hi_wall(false);
703  bool near_y_lo_wall(false); bool near_y_hi_wall(false);
704  bool mask_x_found(false); bool mask_y_found(false);
705 
706  // Near x-wall
707  if ((i-vbx_lo.x) < width) {
708  near_x_lo_wall = true;
709  ii = i-vbx_lo.x;
710  if (mask_arr(vbx_lo.x,j,k) == 2) mask_x_found = true;
711  } else if ((vbx_hi.x-i) < width) {
712  near_x_hi_wall = true;
713  ii = vbx_hi.x-i;
714  if (mask_arr(vbx_hi.x,j,k) == 2) mask_x_found = true;
715  }
716 
717  // Near y-wall
718  if ((j-vbx_lo.y) < width) {
719  near_y_lo_wall = true;
720  jj = j-vbx_lo.y;
721  if (mask_arr(i,vbx_lo.y,k) == 2) mask_y_found = true;
722  } else if ((vbx_hi.y-j) < width) {
723  near_y_hi_wall = true;
724  jj = vbx_hi.y-j;
725  if (mask_arr(i,vbx_hi.y,k) == 2) mask_y_found = true;
726  }
727 
728  // Found a nearby masked cell (valid n_ind)
729  if (mask_x_found && mask_y_found) {
730  n_ind = std::min(ii,jj) + 1.0;
731  } else if (mask_x_found) {
732  n_ind = ii + 1.0;
733  } else if (mask_y_found) {
734  n_ind = jj + 1.0;
735  // Pesky corner cell
736  } else {
737  if (near_x_lo_wall || near_x_hi_wall) {
738  Real dj_min{width-1.0};
739  int j_lb = std::max(vbx_lo.y,j-width);
740  int j_ub = std::min(vbx_hi.y,j+width);
741  int li = (near_x_lo_wall) ? vbx_lo.x : vbx_hi.x;
742  for (int lj(j_lb); lj<=j_ub; ++lj) {
743  if (mask_arr(li,lj,k) == 2) {
744  mask_y_found = true;
745  dj_min = std::min(dj_min,(Real) std::abs(lj-j));
746  }
747  }
748  if (mask_y_found) {
749  Real mag = sqrt( Real(dj_min*dj_min + ii*ii) );
750  n_ind = std::min(mag,width-1.0) + 1.0;
751  } else {
752  Abort("Mask not found near x wall!");
753  }
754  } else if (near_y_lo_wall || near_y_hi_wall) {
755  Real di_min{width-1.0};
756  int i_lb = std::max(vbx_lo.x,i-width);
757  int i_ub = std::min(vbx_hi.x,i+width);
758  int lj = (near_y_lo_wall) ? vbx_lo.y : vbx_hi.y;
759  for (int li(i_lb); li<=i_ub; ++li) {
760  if (mask_arr(li,lj,k) == 2) {
761  mask_x_found = true;
762  di_min = std::min(di_min,(Real) std::abs(li-i));
763  }
764  }
765  if (mask_x_found) {
766  Real mag = sqrt( Real(di_min*di_min + jj*jj) );
767  n_ind = std::min(mag,width-1.0) + 1.0;
768  } else {
769  Abort("Mask not found near y wall!");
770  }
771  } else {
772  Abort("Relaxation cell must be near a wall!");
773  }
774  }
775 
776  Real Factor = (num - n_ind)/denom;
777  Real d = data_arr(i ,j ,k ,n+icomp) + delta_t*rhs_arr(i , j , k ,n+icomp);
778  Real d_ip1 = data_arr(i+1,j ,k ,n+icomp) + delta_t*rhs_arr(i+1, j , k ,n+icomp);
779  Real d_im1 = data_arr(i-1,j ,k ,n+icomp) + delta_t*rhs_arr(i-1, j , k ,n+icomp);
780  Real d_jp1 = data_arr(i ,j+1,k ,n+icomp) + delta_t*rhs_arr(i , j+1, k ,n+icomp);
781  Real d_jm1 = data_arr(i ,j-1,k ,n+icomp) + delta_t*rhs_arr(i , j-1, k ,n+icomp);
782  Real delta = fine_arr(i ,j ,k,n) - d;
783  Real delta_xp = fine_arr(i+1,j ,k,n) - d_ip1;
784  Real delta_xm = fine_arr(i-1,j ,k,n) - d_im1;
785  Real delta_yp = fine_arr(i ,j+1,k,n) - d_jp1;
786  Real delta_ym = fine_arr(i ,j-1,k,n) - d_jm1;
787  Real Laplacian = delta_xp + delta_xm + delta_yp + delta_ym - 4.0*delta;
788  rhs_arr(i,j,k,n) += (F1*delta - F2*Laplacian) * Factor;
789  }
790  });
791  } // mfi
792  } // ivar_idx
793 }
num
@ num
Definition: ERF_DataStruct.H:24
ParallelFor
ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { const auto prob_lo=geomdata.ProbLo();const auto dx=geomdata.CellSize();const Real x=(prob_lo[0]+(i+0.5) *dx[0])/mf_m(i, j, 0);const Real z=z_cc(i, j, k);Real L=std::sqrt(std::pow((x - x_c)/x_r, 2)+std::pow((z - z_c)/z_r, 2));if(L<=1.0) { Real dT=T_pert *(std::cos(PI *L)+1.0)/2.0;Real Tbar_hse=p_hse(i, j, k)/(R_d *r_hse(i, j, k));Real theta_perturbed=(Tbar_hse+dT) *std::pow(p_0/p_hse(i, j, k), rdOcp);Real theta_0=(Tbar_hse) *std::pow(p_0/p_hse(i, j, k), rdOcp);if(const_rho) { state_pert(i, j, k, RhoTheta_comp)=r_hse(i, j, k) *(theta_perturbed - theta_0);} else { state_pert(i, j, k, Rho_comp)=getRhoThetagivenP(p_hse(i, j, k))/theta_perturbed - r_hse(i, j, k);} } })
for
for(int i=0;i< m_num_species;i++)
Definition: ERF_InitCustomPert_MultiSpeciesBubble.H:48
void_bc
PhysBCFunctNoOp void_bc
Definition: ERF_InteriorGhostCells.cpp:5
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:19
rho_arr
auto rho_arr
Definition: ERF_UpdateWSubsidence_SineMassFlux.H:3
ERFFillPatcher::GetMask
amrex::iMultiFab * GetMask()
Definition: ERF_FillPatcher.H:43
ERFFillPatcher::FillRelax
void FillRelax(amrex::MultiFab &mf, amrex::Real time, BC &cbc, amrex::Vector< amrex::BCRec > const &bcs)
Definition: ERF_FillPatcher.H:105
ERFFillPatcher::GetSetMaskVal
int GetSetMaskVal()
Definition: ERF_FillPatcher.H:39
ERFFillPatcher::GetRelaxMaskVal
int GetRelaxMaskVal()
Definition: ERF_FillPatcher.H:41
IntVars::NumTypes
@ NumTypes
Definition: ERF_IndexDefines.H:178
IntVars::ymom
@ ymom
Definition: ERF_IndexDefines.H:176
IntVars::cons
@ cons
Definition: ERF_IndexDefines.H:174
IntVars::zmom
@ zmom
Definition: ERF_IndexDefines.H:177
IntVars::xmom
@ xmom
Definition: ERF_IndexDefines.H:175
Here is the call graph for this function:

◆ realbdy_compute_interior_ghost_rhs()

void realbdy_compute_interior_ghost_rhs ( const Real &  time,
const Real &  delta_t,
const Real &  start_bdy_time,
const Real &  final_bdy_time,
const Real &  bdy_time_interval,
const Real &  nudge_factor,
int  width,
bool  do_upwind,
const Geometry &  geom,
Vector< MultiFab > &  S_rhs,
Vector< MultiFab > &  S_cur_data,
Vector< Vector< FArrayBox >> &  bdy_data_xlo,
Vector< Vector< FArrayBox >> &  bdy_data_xhi,
Vector< Vector< FArrayBox >> &  bdy_data_ylo,
Vector< Vector< FArrayBox >> &  bdy_data_yhi,
std::unique_ptr< ReadBndryPlanes > &  m_r2d 
)

Compute the RHS in the relaxation zone

Parameters
[in]timecurrent (total) time
[in]delta_ttimestep
[in]start_bdy_timefull time of the first time slice of boundary data
[in]final_bdy_timefull time of the last time slice of boundary data
[in]bdy_time_intervaltime interval between boundary condition time stamps
[in]widthnumber of cells in (relaxation+specified) zone
[in]set_widthnumber of cells in (specified) zone
[in]geomcontainer for geometric information
[out]S_rhsRHS to be computed here
[in]S_datacurrent value of the solution
[in]bdy_data_xloboundary data on interior of low x-face
[in]bdy_data_xhiboundary data on interior of high x-face
[in]bdy_data_yloboundary data on interior of low y-face
[in]bdy_data_yhiboundary data on interior of high y-face
112 {
113  BL_PROFILE_REGION("realbdy_compute_interior_ghost_RHS()");
114 
115  // HACK HACK HACK
116  // Get bndry data
117  Vector<int> ind_map2 = {BCVars::xvel_bc, BCVars::yvel_bc, BCVars::RhoTheta_bc_comp};
118  Array4<Real> bdatxlo, bdatxhi, bdatylo, bdatyhi;
119  if (m_r2d) {
120  Vector<std::unique_ptr<PlaneVector>>& bndry_data = m_r2d->interp_in_time(time);
121  bdatxlo = (*bndry_data[0])[0].array();
122  bdatylo = (*bndry_data[1])[0].array();
123  bdatxhi = (*bndry_data[3])[0].array();
124  bdatyhi = (*bndry_data[4])[0].array();
125  }
126 
127  //
128  // Note that time (= start_time+old_stage_time) is measured as total time
129  // start_bdy_time and final_bdy_time are also measured as total time
130  //
131 
132  // Relaxation constants
133  Real F1 = 1./(nudge_factor*delta_t);
134 
135  // Time interpolation
136  Real dT = bdy_time_interval;
137 
138  int n_time = static_cast<int>( (time-start_bdy_time) / dT);
139  int n_time_p1 = n_time + 1;
140  Real alpha = ((time-start_bdy_time) - n_time * dT) / dT;
141 
142  // Do not over run the last bdy file
143  if (time >= final_bdy_time) {
144  n_time = static_cast<int>( (final_bdy_time - start_bdy_time)/ dT);
145  n_time_p1 = n_time;
146  alpha = 0.0;
147  }
148 
149  AMREX_ALWAYS_ASSERT( alpha >= 0. && alpha <= 1.0);
150  Real oma = 1.0 - alpha;
151 
152  /*
153  // UNIT TEST DEBUG
154  oma = 1.0; alpha = 0.0;
155  */
156 
157  // Temporary FABs for storage (owned/filled on all ranks)
158  FArrayBox U_xlo, U_xhi, U_ylo, U_yhi;
159  FArrayBox V_xlo, V_xhi, V_ylo, V_yhi;
160  FArrayBox T_xlo, T_xhi, T_ylo, T_yhi;
161 
162  // Variable index map (WRFBdyVars -> Vars)
163  Vector<int> var_map = {Vars::xvel, Vars::yvel, Vars::cons, Vars::cons};
164  Vector<int> ivar_map = {IntVars::xmom, IntVars::ymom, IntVars::cons, IntVars::cons};
165 
166  // Variable icomp map
167  Vector<int> comp_map = {0, 0, RhoTheta_comp};
168 
169  // Indices
170  int ivarU = RealBdyVars::U;
171  int ivarV = RealBdyVars::V;
172  int ivarT = RealBdyVars::T;
173  int BdyEnd = RealBdyVars::NumTypes-1;
174 
175 
176  // NOTE: The sizing of the temporary BDY FABS is
177  // GLOBAL and occurs over the entire BDY region.
178 
179  // Size the FABs
180  //==========================================================
181  for (int ivar(ivarU); ivar < BdyEnd; ivar++) {
182  int ivar_idx = var_map[ivar];
183  Box domain = geom.Domain();
184  auto ixtype = S_cur_data[ivar_idx].boxArray().ixType();
185  domain.convert(ixtype);
186 
187  // NOTE: Ghost cells needed for idx type mismatch between mask and data (do_upwind)
188  IntVect ng_vect(0);
189  //IntVect ng_vect(1,1,0);
190  Box gdom(domain); gdom.grow(ng_vect);
191  Box bx_xlo, bx_xhi, bx_ylo, bx_yhi;
192  realbdy_interior_bxs_xy(gdom, domain, width,
193  bx_xlo, bx_xhi,
194  bx_ylo, bx_yhi,
195  ng_vect, true);
196 
197  // Size the FABs
198  if (ivar == ivarU) {
199  U_xlo.resize(bx_xlo,1,The_Async_Arena()); U_xhi.resize(bx_xhi,1,The_Async_Arena());
200  U_ylo.resize(bx_ylo,1,The_Async_Arena()); U_yhi.resize(bx_yhi,1,The_Async_Arena());
201  } else if (ivar == ivarV) {
202  V_xlo.resize(bx_xlo,1,The_Async_Arena()); V_xhi.resize(bx_xhi,1,The_Async_Arena());
203  V_ylo.resize(bx_ylo,1,The_Async_Arena()); V_yhi.resize(bx_yhi,1,The_Async_Arena());
204  } else if (ivar == ivarT){
205  T_xlo.resize(bx_xlo,1,The_Async_Arena()); T_xhi.resize(bx_xhi,1,The_Async_Arena());
206  T_ylo.resize(bx_ylo,1,The_Async_Arena()); T_yhi.resize(bx_yhi,1,The_Async_Arena());
207  } else {
208  continue;
209  }
210  } // ivar
211 
212 
213  // NOTE: These operations use the BDY FABS and RHO. The
214  // use of RHO to go from PRIM -> CONS requires that
215  // these operations be LOCAL. So we have allocated
216  // enough space to do global operations (1 rank) but
217  // will fill a subset of that data that the rank owns.
218 
219  // Populate FABs from bdy interpolation (primitive vars)
220  //==========================================================
221  for (int ivar(ivarU); ivar < BdyEnd; ivar++) {
222  int ivar_idx = var_map[ivar];
223  Box domain = geom.Domain();
224  auto ixtype = S_cur_data[ivar_idx].boxArray().ixType();
225  domain.convert(ixtype);
226  const auto& dom_lo = lbound(domain);
227  const auto& dom_hi = ubound(domain);
228 
229  // BndryReg idx and limiting
230  int bdy_comp = ind_map2[ivar];
231  const auto& dom_cc_lo = lbound(geom.Domain());
232  const auto& dom_cc_hi = ubound(geom.Domain());
233 
234 #ifdef _OPENMP
235 #pragma omp parallel if (Gpu::notInLaunchRegion())
236 #endif
237  for (MFIter mfi(S_cur_data[ivar_idx],TilingIfNotGPU()); mfi.isValid(); ++mfi) {
238  // NOTE: Ghost cells needed for idx type mismatch between mask and data (do_upwind)
239  IntVect ng_vect(0);
240  //IntVect ng_vect(1,1,0);
241  Box gtbx = grow(mfi.tilebox(ixtype.toIntVect()),ng_vect);
242  Box tbx_xlo, tbx_xhi, tbx_ylo, tbx_yhi;
243  realbdy_interior_bxs_xy(gtbx, domain, width,
244  tbx_xlo, tbx_xhi,
245  tbx_ylo, tbx_yhi,
246  ng_vect, true);
247 
248  Array4<Real> arr_xlo; Array4<Real> arr_xhi;
249  Array4<Real> arr_ylo; Array4<Real> arr_yhi;
250  if (ivar == ivarU) {
251  arr_xlo = U_xlo.array(); arr_xhi = U_xhi.array();
252  arr_ylo = U_ylo.array(); arr_yhi = U_yhi.array();
253  } else if (ivar == ivarV) {
254  arr_xlo = V_xlo.array(); arr_xhi = V_xhi.array();
255  arr_ylo = V_ylo.array(); arr_yhi = V_yhi.array();
256  } else if (ivar == ivarT){
257  arr_xlo = T_xlo.array(); arr_xhi = T_xhi.array();
258  arr_ylo = T_ylo.array(); arr_yhi = T_yhi.array();
259  } else {
260  continue;
261  }
262 
263  // Boundary data at fixed time intervals
264  const auto& bdatxlo_n = bdy_data_xlo[n_time ][ivar].const_array();
265  const auto& bdatxlo_np1 = bdy_data_xlo[n_time_p1][ivar].const_array();
266  const auto& bdatxhi_n = bdy_data_xhi[n_time ][ivar].const_array();
267  const auto& bdatxhi_np1 = bdy_data_xhi[n_time_p1][ivar].const_array();
268  const auto& bdatylo_n = bdy_data_ylo[n_time ][ivar].const_array();
269  const auto& bdatylo_np1 = bdy_data_ylo[n_time_p1][ivar].const_array();
270  const auto& bdatyhi_n = bdy_data_yhi[n_time ][ivar].const_array();
271  const auto& bdatyhi_np1 = bdy_data_yhi[n_time_p1][ivar].const_array();
272 
273  // Current density to convert to conserved vars
274  Array4<Real> r_arr = S_cur_data[IntVars::cons].array(mfi);
275 
276  // Limiting offset
277  int offset = width - 1;
278 
279  // Populate with interpolation (protect from ghost cells)
280  ParallelFor(tbx_xlo, tbx_xhi,
281  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
282  {
283  int ii = std::max(i , dom_lo.x);
284  ii = std::min(ii, dom_lo.x+offset);
285  int jj = std::max(j , dom_lo.y);
286  jj = std::min(jj, dom_hi.y);
287 
288  Real rho_interp;
289  if (ivar==ivarU) {
290  rho_interp = 0.5 * ( r_arr(i-1,j ,k) + r_arr(i,j,k) );
291  } else if (ivar==ivarV) {
292  rho_interp = 0.5 * ( r_arr(i ,j-1,k) + r_arr(i,j,k) );
293  } else {
294  rho_interp = r_arr(i,j,k);
295  }
296 
297  if (bdatxlo) {
298  int ii2 = std::min(std::max(i , dom_cc_lo.x), dom_cc_hi.x);
299  int jj2 = std::min(std::max(j , dom_cc_lo.y), dom_cc_hi.y);
300  arr_xlo(i,j,k) = rho_interp * bdatxlo(ii2,jj2,k,bdy_comp);
301  } else {
302  arr_xlo(i,j,k) = rho_interp * ( oma * bdatxlo_n (ii,jj,k,0)
303  + alpha * bdatxlo_np1(ii,jj,k,0) );
304  }
305  },
306  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
307  {
308  int ii = std::max(i , dom_hi.x-offset);
309  ii = std::min(ii, dom_hi.x);
310  int jj = std::max(j , dom_lo.y);
311  jj = std::min(jj, dom_hi.y);
312 
313  Real rho_interp;
314  if (ivar==ivarU) {
315  rho_interp = 0.5 * ( r_arr(i-1,j ,k) + r_arr(i,j,k) );
316  } else if (ivar==ivarV) {
317  rho_interp = 0.5 * ( r_arr(i ,j-1,k) + r_arr(i,j,k) );
318  } else {
319  rho_interp = r_arr(i,j,k);
320  }
321 
322  if (bdatxhi) {
323  int ii2 = std::min(std::max(i , dom_cc_lo.x), dom_cc_hi.x);
324  int jj2 = std::min(std::max(j , dom_cc_lo.y), dom_cc_hi.y);
325  arr_xhi(i,j,k) = rho_interp * bdatxhi(ii2,jj2,k,bdy_comp);
326  } else {
327  arr_xhi(i,j,k) = rho_interp * ( oma * bdatxhi_n (ii,jj,k,0)
328  + alpha * bdatxhi_np1(ii,jj,k,0) );
329  }
330  });
331 
332  ParallelFor(tbx_ylo, tbx_yhi,
333  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
334  {
335  int ii = std::max(i , dom_lo.x);
336  ii = std::min(ii, dom_hi.x);
337  int jj = std::max(j , dom_lo.y);
338  jj = std::min(jj, dom_lo.y+offset);
339 
340  Real rho_interp;
341  if (ivar==ivarU) {
342  rho_interp = 0.5 * ( r_arr(i-1,j ,k) + r_arr(i,j,k) );
343  } else if (ivar==ivarV) {
344  rho_interp = 0.5 * ( r_arr(i ,j-1,k) + r_arr(i,j,k) );
345  } else {
346  rho_interp = r_arr(i,j,k);
347  }
348 
349  if (bdatylo) {
350  int ii2 = std::min(std::max(i , dom_cc_lo.x), dom_cc_hi.x);
351  int jj2 = std::min(std::max(j , dom_cc_lo.y), dom_cc_hi.y);
352  arr_ylo(i,j,k) = rho_interp * bdatylo(ii2,jj2,k,bdy_comp);
353  } else {
354  arr_ylo(i,j,k) = rho_interp * ( oma * bdatylo_n (ii,jj,k,0)
355  + alpha * bdatylo_np1(ii,jj,k,0) );
356  }
357  },
358  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
359  {
360  int ii = std::max(i , dom_lo.x);
361  ii = std::min(ii, dom_hi.x);
362  int jj = std::max(j , dom_hi.y-offset);
363  jj = std::min(jj, dom_hi.y);
364 
365  Real rho_interp;
366  if (ivar==ivarU) {
367  rho_interp = 0.5 * ( r_arr(i-1,j ,k) + r_arr(i,j,k) );
368  } else if (ivar==ivarV) {
369  rho_interp = 0.5 * ( r_arr(i ,j-1,k) + r_arr(i,j,k) );
370  } else {
371  rho_interp = r_arr(i,j,k);
372  }
373 
374  if (bdatyhi) {
375  int ii2 = std::min(std::max(i , dom_cc_lo.x), dom_cc_hi.x);
376  int jj2 = std::min(std::max(j , dom_cc_lo.y), dom_cc_hi.y);
377  arr_yhi(i,j,k) = rho_interp * bdatyhi(ii2,jj2,k,bdy_comp);
378  } else {
379  arr_yhi(i,j,k) = rho_interp * ( oma * bdatyhi_n (ii,jj,k,0)
380  + alpha * bdatyhi_np1(ii,jj,k,0) );
381  }
382  });
383  } // mfi
384  } // ivar
385 
386 
387  // Compute RHS in relaxation region
388  //==========================================================
389  auto dx = geom.CellSizeArray();
390  auto ProbLo = geom.ProbLoArray();
391  auto ProbHi = geom.ProbHiArray();
392  for (int ivar(ivarU); ivar < BdyEnd; ivar++) {
393  int ivar_idx = ivar_map[ivar];
394  int icomp = comp_map[ivar];
395 
396  Box domain = geom.Domain();
397  domain.convert(S_cur_data[ivar_idx].boxArray().ixType());
398  IntVect ng_vect(0);
399 
400 #ifdef _OPENMP
401 #pragma omp parallel if (Gpu::notInLaunchRegion())
402 #endif
403  for (MFIter mfi(S_cur_data[ivar_idx],TilingIfNotGPU()); mfi.isValid(); ++mfi) {
404  Box tbx = mfi.tilebox();
405  Box tbx_xlo, tbx_xhi, tbx_ylo, tbx_yhi;
406  realbdy_interior_bxs_xy(tbx, domain, width,
407  tbx_xlo, tbx_xhi,
408  tbx_ylo, tbx_yhi,
409  ng_vect);
410 
411  Array4<Real> rhs_arr; Array4<Real> data_arr;
412  Array4<Real> arr_xlo; Array4<Real> arr_xhi;
413  Array4<Real> arr_ylo; Array4<Real> arr_yhi;
414  if (ivar == ivarU) {
415  arr_xlo = U_xlo.array(); arr_xhi = U_xhi.array();
416  arr_ylo = U_ylo.array(); arr_yhi = U_yhi.array();
417  rhs_arr = S_rhs[IntVars::xmom].array(mfi);
418  data_arr = S_cur_data[IntVars::xmom].array(mfi);
419  } else if (ivar == ivarV) {
420  arr_xlo = V_xlo.array(); arr_xhi = V_xhi.array();
421  arr_ylo = V_ylo.array(); arr_yhi = V_yhi.array();
422  rhs_arr = S_rhs[IntVars::ymom].array(mfi);
423  data_arr = S_cur_data[IntVars::ymom].array(mfi);
424  } else if (ivar == ivarT){
425  arr_xlo = T_xlo.array(); arr_xhi = T_xhi.array();
426  arr_ylo = T_ylo.array(); arr_yhi = T_yhi.array();
427  rhs_arr = S_rhs[IntVars::cons].array(mfi);
428  data_arr = S_cur_data[IntVars::cons].array(mfi);
429  } else {
430  continue;
431  }
432 
433  Array4<Real> u_xlo = U_xlo.array(); Array4<Real> u_xhi = U_xhi.array();
434  Array4<Real> v_xlo = V_xlo.array(); Array4<Real> v_xhi = V_xhi.array();
435  Array4<Real> v_ylo = V_ylo.array(); Array4<Real> v_yhi = V_yhi.array();
436 
437  realbdy_compute_relaxation(icomp, 1,
438  width, dx, ProbLo, ProbHi, F1, geom.Domain(),
439  tbx_xlo , tbx_xhi , tbx_ylo , tbx_yhi ,
440  arr_xlo , arr_xhi , arr_ylo , arr_yhi ,
441  u_xlo, u_xhi, v_xlo, v_xhi, v_ylo, v_yhi,
442  data_arr, rhs_arr, do_upwind);
443 
444  /*
445  // UNIT TEST DEBUG
446  realbdy_interior_bxs_xy(tbx, domain, width,
447  tbx_xlo, tbx_xhi,
448  tbx_ylo, tbx_yhi);
449  ParallelFor(tbx_xlo, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
450  {
451  if (std::fabs(arr_xlo(i,j,k) - data_arr(i,j,k,icomp)) > 0.01) {
452  Print() << "ERROR XLO: " << ivar << ' ' << icomp << ' ' << IntVect(i,j,k) << "\n";
453  Print() << "DATA: " << data_arr(i,j,k,icomp) << ' ' << arr_xlo(i,j,k) << "\n";
454  exit(0);
455  }
456  });
457  ParallelFor(tbx_xhi, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
458  {
459  if (std::fabs(arr_xhi(i,j,k) - data_arr(i,j,k,icomp)) > 0.01) {
460  Print() << "ERROR XHI: " << ivar << ' ' << icomp << ' ' << IntVect(i,j,k) << "\n";
461  Print() << "DATA: " << data_arr(i,j,k,icomp) << ' ' << arr_xhi(i,j,k) << "\n";
462  exit(0);
463  }
464  });
465  ParallelFor(tbx_ylo, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
466  {
467  if (std::fabs(arr_ylo(i,j,k) - data_arr(i,j,k,icomp)) > 0.01) {
468  Print() << "ERROR YLO: " << ivar << ' ' << icomp << ' ' << IntVect(i,j,k) << "\n";
469  Print() << "DATA: " << data_arr(i,j,k,icomp) << ' ' << arr_ylo(i,j,k) << "\n";
470  exit(0);
471  }
472  });
473  ParallelFor(tbx_yhi, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
474  {
475  if (std::fabs(arr_yhi(i,j,k)-data_arr(i,j,k,icomp)) > 0.01) {
476  Print() << "ERROR YHI: " << ivar << ' ' << icomp << ' ' << IntVect(i,j,k) << "\n";
477  Print() << "DATA: " << data_arr(i,j,k,icomp) << ' ' << arr_yhi(i,j,k) << "\n";
478  exit(0);
479  }
480  });
481  */
482  } // mfi
483  } // ivar
484  //ParallelDescriptor::Barrier();
485  //exit(0);
486 }
RhoTheta_comp
#define RhoTheta_comp
Definition: ERF_IndexDefines.H:37
dx
const Real dx
Definition: ERF_InitCustomPert_ABL.H:23
AMREX_ALWAYS_ASSERT
AMREX_ALWAYS_ASSERT(bx.length()[2]==khi+1)
alpha
amrex::Real alpha
Definition: ERF_InitCustomPert_IsentropicVortex.H:8
realbdy_interior_bxs_xy
void realbdy_interior_bxs_xy(const Box &bx, const Box &domain, const int &width, Box &bx_xlo, Box &bx_xhi, Box &bx_ylo, Box &bx_yhi, const IntVect &ng_vect, const bool get_int_ng)
Definition: ERF_InteriorGhostCells.cpp:23
offset
AMREX_FORCE_INLINE IntVect offset(const int face_dir, const int normal)
Definition: ERF_ReadBndryPlanes.cpp:28
dom_hi
const auto & dom_hi
Definition: ERF_SetupVertDiff.H:2
dom_lo
const auto & dom_lo
Definition: ERF_SetupVertDiff.H:1
realbdy_compute_relaxation
AMREX_GPU_HOST AMREX_FORCE_INLINE void realbdy_compute_relaxation(const int &icomp, const int &num_var, const int &width, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &dx, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &ProbLo, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &ProbHi, const amrex::Real &F1, const amrex::Box &domain_cc, const amrex::Box &bx_xlo, const amrex::Box &bx_xhi, const amrex::Box &bx_ylo, const amrex::Box &bx_yhi, const amrex::Array4< const amrex::Real > &arr_xlo, const amrex::Array4< const amrex::Real > &arr_xhi, const amrex::Array4< const amrex::Real > &arr_ylo, const amrex::Array4< const amrex::Real > &arr_yhi, const amrex::Array4< const amrex::Real > &u_xlo, const amrex::Array4< const amrex::Real > &u_xhi, const amrex::Array4< const amrex::Real > &v_xlo, const amrex::Array4< const amrex::Real > &v_xhi, const amrex::Array4< const amrex::Real > &v_ylo, const amrex::Array4< const amrex::Real > &v_yhi, const amrex::Array4< const amrex::Real > &data_arr, const amrex::Array4< amrex::Real > &rhs_arr, const bool &do_upwind)
Definition: ERF_Utils.H:182
BCVars::RhoTheta_bc_comp
@ RhoTheta_bc_comp
Definition: ERF_IndexDefines.H:78
BCVars::yvel_bc
@ yvel_bc
Definition: ERF_IndexDefines.H:88
BCVars::xvel_bc
@ xvel_bc
Definition: ERF_IndexDefines.H:87
RealBdyVars::U
@ U
Definition: ERF_IndexDefines.H:108
RealBdyVars::NumTypes
@ NumTypes
Definition: ERF_IndexDefines.H:112
RealBdyVars::T
@ T
Definition: ERF_IndexDefines.H:110
RealBdyVars::V
@ V
Definition: ERF_IndexDefines.H:109
Vars::xvel
@ xvel
Definition: ERF_IndexDefines.H:157
Vars::cons
@ cons
Definition: ERF_IndexDefines.H:156
Vars::yvel
@ yvel
Definition: ERF_IndexDefines.H:158
Here is the call graph for this function:

◆ realbdy_interior_bxs_xy()

void realbdy_interior_bxs_xy ( const Box &  bx,
const Box &  domain,
const int &  width,
Box &  bx_xlo,
Box &  bx_xhi,
Box &  bx_ylo,
Box &  bx_yhi,
const IntVect &  ng_vect,
const bool  get_int_ng 
)

Get the boxes for looping over interior/exterior ghost cells for use by fillpatch, erf_slow_rhs_pre, and erf_slow_rhs_post.

Parameters
[in]bxbox to intersect with 4 halo regions
[in]domainbox of the whole domain
[in]widthnumber of cells in (relaxation+specified) zone
[in]set_widthnumber of cells in (specified) zone
[out]bx_xlohalo box at x_lo boundary
[out]bx_xhihalo box at x_hi boundary
[out]bx_ylohalo box at y_lo boundary
[out]bx_yhihalo box at y_hi boundary
[in]ng_vectnumber of ghost cells in each direction
[in]get_int_ngflag to get ghost cells inside the domain
32 {
33  AMREX_ALWAYS_ASSERT(bx.ixType() == domain.ixType());
34 
35  //==================================================================
36  // NOTE: X-face boxes take ownership of the overlapping region.
37  // With exterior ghost cells (ng_vect != 0), the x-face
38  // boxes will have exterior ghost cells in both x & y.
39  //==================================================================
40 
41  // Domain bounds without ghost cells
42  const auto& dom_lo = lbound(domain);
43  const auto& dom_hi = ubound(domain);
44 
45  // Four boxes matching the domain
46  Box gdom_xlo(domain); Box gdom_xhi(domain);
47  Box gdom_ylo(domain); Box gdom_yhi(domain);
48 
49  // Trim the boxes to only include internal ghost cells
50  gdom_xlo.setBig(0,dom_lo.x+width-1); gdom_xhi.setSmall(0,dom_hi.x-width+1);
51  gdom_ylo.setBig(1,dom_lo.y+width-1); gdom_yhi.setSmall(1,dom_hi.y-width+1);
52 
53  // Remove overlapping corners from y-face boxes
54  gdom_ylo.setSmall(0,gdom_xlo.bigEnd(0)+1); gdom_ylo.setBig(0,gdom_xhi.smallEnd(0)-1);
55  gdom_yhi.setSmall(0,gdom_xlo.bigEnd(0)+1); gdom_yhi.setBig(0,gdom_xhi.smallEnd(0)-1);
56 
57  // Grow boxes to get external ghost cells only
58  gdom_xlo.growLo(0,ng_vect[0]); gdom_xhi.growHi(0,ng_vect[0]);
59  gdom_xlo.grow (1,ng_vect[1]); gdom_xhi.grow (1,ng_vect[1]);
60  gdom_ylo.growLo(1,ng_vect[1]); gdom_yhi.growHi(1,ng_vect[1]);
61 
62  // Grow boxes to get internal ghost cells
63  if (get_int_ng) {
64  gdom_xlo.growHi(0,ng_vect[0]); gdom_xhi.growLo(0,ng_vect[0]);
65  gdom_ylo.grow (0,ng_vect[0]); gdom_yhi.grow (0,ng_vect[0]);
66  gdom_ylo.growHi(1,ng_vect[1]); gdom_yhi.growLo(1,ng_vect[1]);
67  }
68 
69  // Populate everything
70  bx_xlo = (bx & gdom_xlo);
71  bx_xhi = (bx & gdom_xhi);
72  bx_ylo = (bx & gdom_ylo);
73  bx_yhi = (bx & gdom_yhi);
74 }

Referenced by realbdy_compute_interior_ghost_rhs().

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

Variable Documentation

◆ void_bc

PhysBCFunctNoOp void_bc

Referenced by fine_compute_interior_ghost_rhs().