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_bc_bxs_xy (const Box &bx, const Box &domain, const int &set_width, Box &bx_xlo, Box &bx_xhi, Box &bx_ylo, Box &bx_yhi, const IntVect &ng_vect)
 
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, 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, const Real &c_p, const Real &rdOcp)
 
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
644 {
645  BL_PROFILE_REGION("fine_compute_interior_ghost_RHS()");
646 
647  // Relaxation constants
648  Real F1 = one/(Real(10.)*delta_t);
649  Real F2 = one/(Real(50.)*delta_t);
650 
651  // Vector of MFs to hold data (dm differs w/ fine patch)
652  Vector<MultiFab> fmf_p_v;
653 
654  // Loop over the variables
655  for (int ivar_idx = 0; ivar_idx < IntVars::NumTypes; ++ivar_idx)
656  {
657  // Fine mfs
658  MultiFab& fmf = S_data_f[ivar_idx];
659  MultiFab& rhs = S_rhs_f [ivar_idx];
660 
661  // NOTE: These temporary MFs and copy operations are horrible
662  // for memory usage and efficiency. However, we need to
663  // have access to ghost cells in the cons array to convert
664  // from primitive u/v/w to momentum. Furthermore, the BA
665  // for the fine patches in ERFFillPatcher don't match the
666  // BA for the data/RHS. For this reason, the data is copied
667  // to a vector of MFs (with ghost cells) so the BAs match
668  // the BA of data/RHS and we have access to rho to convert
669  // prim to conserved.
670 
671  // Temp MF on box (distribution map differs w/ fine patch)
672  int num_var = fmf.nComp();
673  fmf_p_v.emplace_back(fmf.boxArray(), fmf.DistributionMap(), num_var, fmf.nGrowVect());
674  MultiFab& fmf_p = fmf_p_v[ivar_idx];
675  MultiFab::Copy(fmf_p,fmf, 0, 0, num_var, fmf.nGrowVect());
676 
677  // Integer mask MF
678  int set_mask_val;
679  int relax_mask_val;
680  iMultiFab* mask;
681 
682  // Fill fine patch on interior halo region
683  //==========================================================
684  if (ivar_idx == IntVars::cons)
685  {
686  FPr_c->FillRelax(fmf_p, time, void_bc, domain_bcs_type);
687  mask = FPr_c->GetMask();
688  set_mask_val = FPr_c->GetSetMaskVal();
689  relax_mask_val = FPr_c->GetRelaxMaskVal();
690  }
691  else if (ivar_idx == IntVars::xmom)
692  {
693  FPr_u->FillRelax(fmf_p, time, void_bc, domain_bcs_type);
694  mask = FPr_u->GetMask();
695  set_mask_val = FPr_u->GetSetMaskVal();
696  relax_mask_val = FPr_u->GetRelaxMaskVal();
697 
698 #ifdef _OPENMP
699 #pragma omp parallel if (Gpu::notInLaunchRegion())
700 #endif
701  for ( MFIter mfi(fmf_p,TilingIfNotGPU()); mfi.isValid(); ++mfi)
702  {
703  Box tbx = mfi.tilebox();
704 
705  const Array4<Real>& prim_arr = fmf_p.array(mfi);
706  const Array4<const Real>& rho_arr = fmf_p_v[0].const_array(mfi);
707  const Array4<const int>& mask_arr = mask->const_array(mfi);
708 
709  ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
710  {
711  if (mask_arr(i,j,k) == relax_mask_val) {
712  Real rho_interp = myhalf * ( rho_arr(i-1,j,k) + rho_arr(i,j,k) );
713  prim_arr(i,j,k) *= rho_interp;
714  }
715  });
716  } // mfi
717  }
718  else if (ivar_idx == IntVars::ymom)
719  {
720  FPr_v->FillRelax(fmf_p, time, void_bc, domain_bcs_type);
721  mask = FPr_v->GetMask();
722  set_mask_val = FPr_v->GetSetMaskVal();
723  relax_mask_val = FPr_v->GetRelaxMaskVal();
724 
725 #ifdef _OPENMP
726 #pragma omp parallel if (Gpu::notInLaunchRegion())
727 #endif
728  for ( MFIter mfi(fmf_p,TilingIfNotGPU()); mfi.isValid(); ++mfi)
729  {
730  Box tbx = mfi.tilebox();
731 
732  const Array4<Real>& prim_arr = fmf_p.array(mfi);
733  const Array4<const Real>& rho_arr = fmf_p_v[0].const_array(mfi);
734  const Array4<const int>& mask_arr = mask->const_array(mfi);
735 
736  ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
737  {
738  if (mask_arr(i,j,k) == relax_mask_val) {
739  Real rho_interp = myhalf * ( rho_arr(i,j-1,k) + rho_arr(i,j,k) );
740  prim_arr(i,j,k) *= rho_interp;
741  }
742  });
743  } // mfi
744  }
745  else if (ivar_idx == IntVars::zmom)
746  {
747  FPr_w->FillRelax(fmf_p, time, void_bc, domain_bcs_type);
748  mask = FPr_w->GetMask();
749  set_mask_val = FPr_w->GetSetMaskVal();
750  relax_mask_val = FPr_w->GetRelaxMaskVal();
751 
752 #ifdef _OPENMP
753 #pragma omp parallel if (Gpu::notInLaunchRegion())
754 #endif
755  for ( MFIter mfi(fmf_p,TilingIfNotGPU()); mfi.isValid(); ++mfi)
756  {
757  Box tbx = mfi.tilebox();
758 
759  const Array4<Real>& prim_arr = fmf_p.array(mfi);
760  const Array4<const Real>& rho_arr = fmf_p_v[0].const_array(mfi);
761  const Array4<const int>& mask_arr = mask->const_array(mfi);
762 
763  ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
764  {
765  if (mask_arr(i,j,k) == relax_mask_val) {
766  Real rho_interp = myhalf * ( rho_arr(i,j,k-1) + rho_arr(i,j,k) );
767  prim_arr(i,j,k) *= rho_interp;
768  }
769  });
770  } // mfi
771  } else {
772  Abort("Dont recognize this variable type in fine_compute_interior_ghost_RHS");
773  }
774 
775 
776  // Zero RHS in set region
777  //==========================================================
778 #ifdef _OPENMP
779 #pragma omp parallel if (Gpu::notInLaunchRegion())
780 #endif
781  for ( MFIter mfi(rhs,TilingIfNotGPU()); mfi.isValid(); ++mfi)
782  {
783  Box tbx = mfi.tilebox();
784  const Array4<Real>& rhs_arr = rhs.array(mfi);
785  const Array4<const int>& mask_arr = mask->const_array(mfi);
786 
787  ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
788  {
789  if (mask_arr(i,j,k) == set_mask_val) {
790  rhs_arr(i,j,k) = zero;
791  }
792  });
793  } // mfi
794 
795  // For Laplacian stencil
796  rhs.FillBoundary(geom.periodicity());
797 
798 
799  // Compute RHS in relaxation region
800  //==========================================================
801 #ifdef _OPENMP
802 #pragma omp parallel if (Gpu::notInLaunchRegion())
803 #endif
804  for ( MFIter mfi(fmf_p,TilingIfNotGPU()); mfi.isValid(); ++mfi)
805  {
806  Box tbx = mfi.tilebox();
807  const Array4<Real>& rhs_arr = rhs.array(mfi);
808  const Array4<const Real>& fine_arr = fmf_p.const_array(mfi);
809  const Array4<const Real>& data_arr = fmf.const_array(mfi);
810  const Array4<const int>& mask_arr = mask->const_array(mfi);
811 
812  Box vbx = mfi.validbox();
813  const auto& vbx_lo = lbound(vbx);
814  const auto& vbx_hi = ubound(vbx);
815 
816  int icomp = 0;
817 
818  int Spec_z = set_width;
819  int Relax_z = width - Spec_z;
820  Real num = Real(Spec_z + Relax_z);
821  Real denom = Real(Relax_z - 1);
822  ParallelFor(tbx, num_var, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
823  {
824  if (mask_arr(i,j,k) == relax_mask_val) {
825 
826  // Indices
827  Real n_ind(-1); // Set to -1 to quiet compiler warning
828  int ii(width-1); int jj(width-1);
829  bool near_x_lo_wall(false); bool near_x_hi_wall(false);
830  bool near_y_lo_wall(false); bool near_y_hi_wall(false);
831  bool mask_x_found(false); bool mask_y_found(false);
832 
833  // Near x-wall
834  if ((i-vbx_lo.x) < width) {
835  near_x_lo_wall = true;
836  ii = i-vbx_lo.x;
837  if (mask_arr(vbx_lo.x,j,k) == 2) mask_x_found = true;
838  } else if ((vbx_hi.x-i) < width) {
839  near_x_hi_wall = true;
840  ii = vbx_hi.x-i;
841  if (mask_arr(vbx_hi.x,j,k) == 2) mask_x_found = true;
842  }
843 
844  // Near y-wall
845  if ((j-vbx_lo.y) < width) {
846  near_y_lo_wall = true;
847  jj = j-vbx_lo.y;
848  if (mask_arr(i,vbx_lo.y,k) == 2) mask_y_found = true;
849  } else if ((vbx_hi.y-j) < width) {
850  near_y_hi_wall = true;
851  jj = vbx_hi.y-j;
852  if (mask_arr(i,vbx_hi.y,k) == 2) mask_y_found = true;
853  }
854 
855  // Found a nearby masked cell (valid n_ind)
856  if (mask_x_found && mask_y_found) {
857  n_ind = std::min(ii,jj) + one;
858  } else if (mask_x_found) {
859  n_ind = ii + one;
860  } else if (mask_y_found) {
861  n_ind = jj + one;
862  // Pesky corner cell
863  } else {
864  if (near_x_lo_wall || near_x_hi_wall) {
865  Real dj_min{width-one};
866  int j_lb = std::max(vbx_lo.y,j-width);
867  int j_ub = std::min(vbx_hi.y,j+width);
868  int li = (near_x_lo_wall) ? vbx_lo.x : vbx_hi.x;
869  for (int lj(j_lb); lj<=j_ub; ++lj) {
870  if (mask_arr(li,lj,k) == 2) {
871  mask_y_found = true;
872  dj_min = std::min(dj_min,(Real) std::abs(lj-j));
873  }
874  }
875  if (mask_y_found) {
876  Real mag = std::sqrt( Real(dj_min*dj_min + ii*ii) );
877  n_ind = std::min(mag,width-one) + one;
878  } else {
879  Abort("Mask not found near x wall!");
880  }
881  } else if (near_y_lo_wall || near_y_hi_wall) {
882  Real di_min{width-one};
883  int i_lb = std::max(vbx_lo.x,i-width);
884  int i_ub = std::min(vbx_hi.x,i+width);
885  int lj = (near_y_lo_wall) ? vbx_lo.y : vbx_hi.y;
886  for (int li(i_lb); li<=i_ub; ++li) {
887  if (mask_arr(li,lj,k) == 2) {
888  mask_x_found = true;
889  di_min = std::min(di_min,(Real) std::abs(li-i));
890  }
891  }
892  if (mask_x_found) {
893  Real mag = std::sqrt( Real(di_min*di_min + jj*jj) );
894  n_ind = std::min(mag,width-one) + one;
895  } else {
896  Abort("Mask not found near y wall!");
897  }
898  } else {
899  Abort("Relaxation cell must be near a wall!");
900  }
901  }
902 
903  Real Factor = (num - n_ind)/denom;
904  Real d = data_arr(i ,j ,k ,n+icomp) + delta_t*rhs_arr(i , j , k ,n+icomp);
905  Real d_ip1 = data_arr(i+1,j ,k ,n+icomp) + delta_t*rhs_arr(i+1, j , k ,n+icomp);
906  Real d_im1 = data_arr(i-1,j ,k ,n+icomp) + delta_t*rhs_arr(i-1, j , k ,n+icomp);
907  Real d_jp1 = data_arr(i ,j+1,k ,n+icomp) + delta_t*rhs_arr(i , j+1, k ,n+icomp);
908  Real d_jm1 = data_arr(i ,j-1,k ,n+icomp) + delta_t*rhs_arr(i , j-1, k ,n+icomp);
909  Real delta = fine_arr(i ,j ,k,n) - d;
910  Real delta_xp = fine_arr(i+1,j ,k,n) - d_ip1;
911  Real delta_xm = fine_arr(i-1,j ,k,n) - d_im1;
912  Real delta_yp = fine_arr(i ,j+1,k,n) - d_jp1;
913  Real delta_ym = fine_arr(i ,j-1,k,n) - d_jm1;
914  Real Laplacian = delta_xp + delta_xm + delta_yp + delta_ym - Real(4.0)*delta;
915  rhs_arr(i,j,k,n) += (F1*delta - F2*Laplacian) * Factor;
916  }
917  });
918  } // mfi
919  } // ivar_idx
920 }
one
constexpr amrex::Real one
Definition: ERF_Constants.H:9
zero
constexpr amrex::Real zero
Definition: ERF_Constants.H:8
myhalf
constexpr amrex::Real myhalf
Definition: ERF_Constants.H:13
num
@ num
Definition: ERF_DataStruct.H:24
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
ParallelFor
ParallelFor(grown_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);})
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:196
IntVars::ymom
@ ymom
Definition: ERF_IndexDefines.H:194
IntVars::cons
@ cons
Definition: ERF_IndexDefines.H:192
IntVars::zmom
@ zmom
Definition: ERF_IndexDefines.H:195
IntVars::xmom
@ xmom
Definition: ERF_IndexDefines.H:193
Here is the call graph for this function:

◆ realbdy_bc_bxs_xy()

void realbdy_bc_bxs_xy ( const Box &  bx,
const Box &  domain,
const int &  set_width,
Box &  bx_xlo,
Box &  bx_xhi,
Box &  bx_ylo,
Box &  bx_yhi,
const IntVect &  ng_vect 
)

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
101 {
102  AMREX_ALWAYS_ASSERT(bx.ixType() == domain.ixType());
103 
104  // Domain bounds without ghost cells
105  const auto& dom_lo = lbound(domain);
106  const auto& dom_hi = ubound(domain);
107 
108  // Four boxes matching the domain
109  Box gdom_xlo(domain); Box gdom_xhi(domain);
110  Box gdom_ylo(domain); Box gdom_yhi(domain);
111 
112  // Get offsets from box index type
113  IntVect iv_type = bx.ixType().toIntVect();
114  int offx = (iv_type[0]==1) ? 0 : -1;
115  int offy = (iv_type[1]==1) ? 0 : -1;
116 
117  // Stagger the boxes based upon index type
118  gdom_xlo += IntVect(offx,0,0); gdom_xhi += IntVect(-offx,0,0);
119  gdom_ylo += IntVect(0,offy,0); gdom_yhi += IntVect(0,-offy,0);
120 
121  // Trim the boxes to only include internal ghost cells
122  gdom_xlo.setBig(0,dom_lo.x+set_width+offx-1); gdom_xhi.setSmall(0,dom_hi.x-set_width-offx+1);
123  gdom_ylo.setBig(1,dom_lo.y+set_width+offy-1); gdom_yhi.setSmall(1,dom_hi.y-set_width-offy+1);
124 
125  // Remove overlapping corners from y-face boxes
126  gdom_ylo.setSmall(0,gdom_xlo.bigEnd(0)+1); gdom_ylo.setBig(0,gdom_xhi.smallEnd(0)-1);
127  gdom_yhi.setSmall(0,gdom_xlo.bigEnd(0)+1); gdom_yhi.setBig(0,gdom_xhi.smallEnd(0)-1);
128 
129  // Grow boxes to get external ghost cells only
130  gdom_xlo.growLo(0,ng_vect[0]+offx); gdom_xhi.growHi(0,ng_vect[0]+offx);
131  gdom_xlo.grow (1,ng_vect[1] ); gdom_xhi.grow (1,ng_vect[1] );
132  gdom_ylo.growLo(1,ng_vect[1]+offy); gdom_yhi.growHi(1,ng_vect[1]+offy);
133 
134  // Populate everything
135  bx_xlo = (bx & gdom_xlo);
136  bx_xhi = (bx & gdom_xhi);
137  bx_ylo = (bx & gdom_ylo);
138  bx_yhi = (bx & gdom_yhi);
139 }
AMREX_ALWAYS_ASSERT
AMREX_ALWAYS_ASSERT(bx.length()[2]==khi+1)
dom_hi
const auto & dom_hi
Definition: ERF_SetupVertDiff.H:2
dom_lo
const auto & dom_lo
Definition: ERF_SetupVertDiff.H:1
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,
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,
const Real &  c_p,
const Real &  rdOcp 
)

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
178 {
179  BL_PROFILE_REGION("realbdy_compute_interior_ghost_RHS()");
180 
181  //
182  // Note that time (= start_time+old_stage_time) is measured as total time
183  // start_bdy_time and final_bdy_time are also measured as total time
184  //
185 
186  // Get bndry data if we have it
187  Vector<int> bnd_map = {BCVars::xvel_bc, BCVars::yvel_bc, BCVars::RhoTheta_bc_comp};
188  Array4<Real> bdatxlo, bdatxhi, bdatylo, bdatyhi;
189  Array4<Real> btenxlo, btenxhi, btenylo, btenyhi;
190  if (m_r2d) {
191  // Index is [plane orientation] and [level]
192  Vector<std::unique_ptr<PlaneVector>>& bndry_data = m_r2d->interp_in_time(time);
193  bdatxlo = (*bndry_data[0])[0].array();
194  bdatylo = (*bndry_data[1])[0].array();
195  bdatxhi = (*bndry_data[3])[0].array();
196  bdatyhi = (*bndry_data[4])[0].array();
197 
198  Vector<std::unique_ptr<PlaneVector>>& bndry_tend = m_r2d->get_tendency(time);
199  btenxlo = (*bndry_tend[0])[0].array();
200  btenylo = (*bndry_tend[1])[0].array();
201  btenxhi = (*bndry_tend[3])[0].array();
202  btenyhi = (*bndry_tend[4])[0].array();
203  }
204 
205  // Relaxation constants
206  Real F1 = one/(nudge_factor*delta_t);
207 
208  // Time interpolation
209  Real dT = bdy_time_interval;
210 
211  int n_time = static_cast<int>( (time-start_bdy_time) / dT);
212  int n_time_p1 = n_time + 1;
213  Real alpha = ((time-start_bdy_time) - n_time * dT) / dT;
214 
215  // Do not over run the last bdy file
216  if (time >= final_bdy_time) {
217  n_time = static_cast<int>( (final_bdy_time - start_bdy_time)/ dT);
218  n_time_p1 = n_time;
219  alpha = zero;
220  }
221 
222  AMREX_ALWAYS_ASSERT( alpha >= zero && alpha <= one);
223  Real oma = one - alpha;
224 
225  // Temporary FABs for storage (owned/filled on all ranks)
226  FArrayBox U_xlo, U_xhi, U_ylo, U_yhi;
227  FArrayBox V_xlo, V_xhi, V_ylo, V_yhi;
228  FArrayBox T_xlo, T_xhi, T_ylo, T_yhi;
229 
230  // Variable index map (WRFBdyVars -> Vars)
231  Vector<int> var_map = {Vars::xvel, Vars::yvel, Vars::cons };
232  Vector<int> ivar_map = {IntVars::xmom, IntVars::ymom, IntVars::cons};
233 
234  // Variable icomp map
235  Vector<int> comp_map = {0, 0, RhoTheta_comp};
236 
237  // Indices
238  int ivarU = RealBdyVars::U;
239  int ivarV = RealBdyVars::V;
240  int ivarT = RealBdyVars::T;
241  int BdyEnd = RealBdyVars::NumTypes-1;
242 
243 
244  // NOTE: The sizing of the temporary BDY FABS is
245  // GLOBAL and occurs over the entire BDY region.
246 
247  // Size the FABs
248  //==========================================================
249  for (int ivar(ivarU); ivar < BdyEnd; ivar++) {
250  int ivar_idx = var_map[ivar];
251  Box domain = geom.Domain();
252  auto ixtype = S_cur_data[ivar_idx].boxArray().ixType();
253  domain.convert(ixtype);
254 
255  // NOTE: Ghost cells needed for idx type mismatch between mask and data (do_upwind)
256  IntVect ng_vect(0);
257  //IntVect ng_vect(1,1,0);
258  Box gdom(domain); gdom.grow(ng_vect);
259  Box bx_xlo, bx_xhi, bx_ylo, bx_yhi;
260  realbdy_interior_bxs_xy(gdom, domain, width,
261  bx_xlo, bx_xhi,
262  bx_ylo, bx_yhi,
263  ng_vect, true);
264 
265  // Size the FABs
266  if (ivar == ivarU) {
267  U_xlo.resize(bx_xlo,1,The_Async_Arena()); U_xhi.resize(bx_xhi,1,The_Async_Arena());
268  U_ylo.resize(bx_ylo,1,The_Async_Arena()); U_yhi.resize(bx_yhi,1,The_Async_Arena());
269  } else if (ivar == ivarV) {
270  V_xlo.resize(bx_xlo,1,The_Async_Arena()); V_xhi.resize(bx_xhi,1,The_Async_Arena());
271  V_ylo.resize(bx_ylo,1,The_Async_Arena()); V_yhi.resize(bx_yhi,1,The_Async_Arena());
272  } else if (ivar == ivarT){
273  T_xlo.resize(bx_xlo,1,The_Async_Arena()); T_xhi.resize(bx_xhi,1,The_Async_Arena());
274  T_ylo.resize(bx_ylo,1,The_Async_Arena()); T_yhi.resize(bx_yhi,1,The_Async_Arena());
275  } else {
276  continue;
277  }
278  } // ivar
279 
280 
281  // NOTE: These operations use the BDY FABS and RHO. The
282  // use of RHO to go from PRIM -> CONS requires that
283  // these operations be LOCAL. So we have allocated
284  // enough space to do global operations (1 rank) but
285  // will fill a subset of that data that the rank owns.
286 
287  // Populate FABs from bdy interpolation (primitive vars)
288  //==========================================================
289  for (int ivar(ivarU); ivar < BdyEnd; ivar++) {
290  int ivar_idx = var_map[ivar];
291  Box domain = geom.Domain();
292  auto ixtype = S_cur_data[ivar_idx].boxArray().ixType();
293  domain.convert(ixtype);
294  const auto& dom_lo = lbound(domain);
295  const auto& dom_hi = ubound(domain);
296 
297  // BndryReg idx and limiting
298  int bdy_comp = bnd_map[ivar];
299  const auto& dom_cc_lo = lbound(geom.Domain());
300  const auto& dom_cc_hi = ubound(geom.Domain());
301 
302 #ifdef _OPENMP
303 #pragma omp parallel if (Gpu::notInLaunchRegion())
304 #endif
305  for (MFIter mfi(S_cur_data[ivar_idx],TilingIfNotGPU()); mfi.isValid(); ++mfi) {
306  // NOTE: Ghost cells needed for idx type mismatch between mask and data (do_upwind)
307  IntVect ng_vect(0);
308  //IntVect ng_vect(1,1,0);
309  Box gtbx = grow(mfi.tilebox(ixtype.toIntVect()),ng_vect);
310  Box tbx_xlo, tbx_xhi, tbx_ylo, tbx_yhi;
311  realbdy_interior_bxs_xy(gtbx, domain, width,
312  tbx_xlo, tbx_xhi,
313  tbx_ylo, tbx_yhi,
314  ng_vect, true);
315 
316  Array4<Real> arr_xlo; Array4<Real> arr_xhi;
317  Array4<Real> arr_ylo; Array4<Real> arr_yhi;
318  if (ivar == ivarU) {
319  arr_xlo = U_xlo.array(); arr_xhi = U_xhi.array();
320  arr_ylo = U_ylo.array(); arr_yhi = U_yhi.array();
321  } else if (ivar == ivarV) {
322  arr_xlo = V_xlo.array(); arr_xhi = V_xhi.array();
323  arr_ylo = V_ylo.array(); arr_yhi = V_yhi.array();
324  } else if (ivar == ivarT){
325  arr_xlo = T_xlo.array(); arr_xhi = T_xhi.array();
326  arr_ylo = T_ylo.array(); arr_yhi = T_yhi.array();
327  } else {
328  continue;
329  }
330 
331  // Boundary data at fixed time intervals
332  const auto& bdatxlo_n = bdy_data_xlo[n_time ][ivar].const_array();
333  const auto& bdatxlo_np1 = bdy_data_xlo[n_time_p1][ivar].const_array();
334  const auto& bdatxhi_n = bdy_data_xhi[n_time ][ivar].const_array();
335  const auto& bdatxhi_np1 = bdy_data_xhi[n_time_p1][ivar].const_array();
336  const auto& bdatylo_n = bdy_data_ylo[n_time ][ivar].const_array();
337  const auto& bdatylo_np1 = bdy_data_ylo[n_time_p1][ivar].const_array();
338  const auto& bdatyhi_n = bdy_data_yhi[n_time ][ivar].const_array();
339  const auto& bdatyhi_np1 = bdy_data_yhi[n_time_p1][ivar].const_array();
340 
341  // Current density to convert to conserved vars
342  Array4<Real> r_arr = S_cur_data[IntVars::cons].array(mfi);
343 
344  // Limiting offset
345  int offset = width - 1;
346 
347  // Populate with interpolation (protect from ghost cells)
348  ParallelFor(tbx_xlo, tbx_xhi,
349  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
350  {
351  int ii = std::max(i , dom_lo.x); ii = std::min(ii, dom_lo.x+offset);
352  int jj = std::max(j , dom_lo.y); jj = std::min(jj, dom_hi.y);
353 
354  Real rho_interp;
355  if (ivar==ivarU) {
356  rho_interp = myhalf * ( r_arr(i-1,j ,k) + r_arr(i,j,k) );
357  } else if (ivar==ivarV) {
358  rho_interp = myhalf * ( r_arr(i ,j-1,k) + r_arr(i,j,k) );
359  } else {
360  rho_interp = r_arr(i,j,k);
361  }
362 
363  if (bdatxlo) {
364  int ii2 = std::min(std::max(i , dom_cc_lo.x), dom_cc_hi.x);
365  int jj2 = std::min(std::max(j , dom_cc_lo.y), dom_cc_hi.y);
366  arr_xlo(i,j,k) = rho_interp * bdatxlo(ii2,jj2,k,bdy_comp);
367  } else {
368  arr_xlo(i,j,k) = rho_interp * ( oma * bdatxlo_n (ii,jj,k,0)
369  + alpha * bdatxlo_np1(ii,jj,k,0) );
370  }
371  },
372  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
373  {
374  int ii = std::max(i , dom_hi.x-offset); ii = std::min(ii, dom_hi.x);
375  int jj = std::max(j , dom_lo.y); jj = std::min(jj, dom_hi.y);
376 
377  Real rho_interp;
378  if (ivar==ivarU) {
379  rho_interp = myhalf * ( r_arr(i-1,j ,k) + r_arr(i,j,k) );
380  } else if (ivar==ivarV) {
381  rho_interp = myhalf * ( r_arr(i ,j-1,k) + r_arr(i,j,k) );
382  } else {
383  rho_interp = r_arr(i,j,k);
384  }
385 
386  if (bdatxhi) {
387  int ii2 = std::min(std::max(i , dom_cc_lo.x), dom_cc_hi.x);
388  int jj2 = std::min(std::max(j , dom_cc_lo.y), dom_cc_hi.y);
389  arr_xhi(i,j,k) = rho_interp * bdatxhi(ii2,jj2,k,bdy_comp);
390  } else {
391  arr_xhi(i,j,k) = rho_interp * ( oma * bdatxhi_n (ii,jj,k,0)
392  + alpha * bdatxhi_np1(ii,jj,k,0) );
393  }
394  });
395 
396  ParallelFor(tbx_ylo, tbx_yhi,
397  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
398  {
399  int ii = std::max(i , dom_lo.x); ii = std::min(ii, dom_hi.x);
400  int jj = std::max(j , dom_lo.y); jj = std::min(jj, dom_lo.y+offset);
401 
402  Real rho_interp;
403  if (ivar==ivarU) {
404  rho_interp = myhalf * ( r_arr(i-1,j ,k) + r_arr(i,j,k) );
405  } else if (ivar==ivarV) {
406  rho_interp = myhalf * ( r_arr(i ,j-1,k) + r_arr(i,j,k) );
407  } else {
408  rho_interp = r_arr(i,j,k);
409  }
410 
411  if (bdatylo) {
412  int ii2 = std::min(std::max(i , dom_cc_lo.x), dom_cc_hi.x);
413  int jj2 = std::min(std::max(j , dom_cc_lo.y), dom_cc_hi.y);
414  arr_ylo(i,j,k) = rho_interp * bdatylo(ii2,jj2,k,bdy_comp);
415  } else {
416  arr_ylo(i,j,k) = rho_interp * ( oma * bdatylo_n (ii,jj,k,0)
417  + alpha * bdatylo_np1(ii,jj,k,0) );
418  }
419  },
420  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
421  {
422  int ii = std::max(i , dom_lo.x); ii = std::min(ii, dom_hi.x);
423  int jj = std::max(j , dom_hi.y-offset); jj = std::min(jj, dom_hi.y);
424 
425  Real rho_interp;
426  if (ivar==ivarU) {
427  rho_interp = myhalf * ( r_arr(i-1,j ,k) + r_arr(i,j,k) );
428  } else if (ivar==ivarV) {
429  rho_interp = myhalf * ( r_arr(i ,j-1,k) + r_arr(i,j,k) );
430  } else {
431  rho_interp = r_arr(i,j,k);
432  }
433 
434  if (bdatyhi) {
435  int ii2 = std::min(std::max(i , dom_cc_lo.x), dom_cc_hi.x);
436  int jj2 = std::min(std::max(j , dom_cc_lo.y), dom_cc_hi.y);
437  arr_yhi(i,j,k) = rho_interp * bdatyhi(ii2,jj2,k,bdy_comp);
438  } else {
439  arr_yhi(i,j,k) = rho_interp * ( oma * bdatyhi_n (ii,jj,k,0)
440  + alpha * bdatyhi_np1(ii,jj,k,0) );
441  }
442  });
443  } // mfi
444  } // ivar
445 
446 
447  // Compute RHS in relaxation region
448  //==========================================================
449  auto dx = geom.CellSizeArray();
450  auto ProbLo = geom.ProbLoArray();
451  auto ProbHi = geom.ProbHiArray();
452 
453  for (int ivar(ivarU); ivar < BdyEnd; ivar++) {
454  int ivar_idx = ivar_map[ivar];
455  int icomp = comp_map[ivar];
456 
457  Box domain = geom.Domain();
458  domain.convert(S_cur_data[ivar_idx].boxArray().ixType());
459  IntVect ng_vect(0);
460 
461 #ifdef _OPENMP
462 #pragma omp parallel if (Gpu::notInLaunchRegion())
463 #endif
464  for (MFIter mfi(S_cur_data[ivar_idx],TilingIfNotGPU()); mfi.isValid(); ++mfi)
465  {
466  Box tbx = mfi.tilebox();
467  Box tbx_xlo, tbx_xhi, tbx_ylo, tbx_yhi;
468  realbdy_interior_bxs_xy(tbx, domain, width,
469  tbx_xlo, tbx_xhi,
470  tbx_ylo, tbx_yhi,
471  ng_vect);
472 
473  Array4<Real> rhs_arr; Array4<Real> data_arr;
474  Array4<Real> arr_xlo; Array4<Real> arr_xhi;
475  Array4<Real> arr_ylo; Array4<Real> arr_yhi;
476  if (ivar == ivarU) {
477  arr_xlo = U_xlo.array(); arr_xhi = U_xhi.array();
478  arr_ylo = U_ylo.array(); arr_yhi = U_yhi.array();
479  rhs_arr = S_rhs[IntVars::xmom].array(mfi);
480  data_arr = S_cur_data[IntVars::xmom].array(mfi);
481  } else if (ivar == ivarV) {
482  arr_xlo = V_xlo.array(); arr_xhi = V_xhi.array();
483  arr_ylo = V_ylo.array(); arr_yhi = V_yhi.array();
484  rhs_arr = S_rhs[IntVars::ymom].array(mfi);
485  data_arr = S_cur_data[IntVars::ymom].array(mfi);
486  } else if (ivar == ivarT){
487  arr_xlo = T_xlo.array(); arr_xhi = T_xhi.array();
488  arr_ylo = T_ylo.array(); arr_yhi = T_yhi.array();
489  rhs_arr = S_rhs[IntVars::cons].array(mfi);
490  data_arr = S_cur_data[IntVars::cons].array(mfi);
491  } else {
492  continue;
493  }
494 
495  realbdy_compute_relaxation(icomp, 1,
496  width, dx, ProbLo, ProbHi, F1,
497  tbx_xlo , tbx_xhi , tbx_ylo , tbx_yhi ,
498  arr_xlo , arr_xhi , arr_ylo , arr_yhi ,
499  data_arr, rhs_arr , c_p, rdOcp);
500  } // mfi
501  } // ivar
502 
503  // Set normal velocity RHS at the boundary
504  //==========================================================
505  Box domain = geom.Domain();
506  Box domainx = convert(domain, IntVect(1,0,0));
507  Box domainy = convert(domain, IntVect(0,1,0));
508 
509  int ilo = domainx.smallEnd(0);
510  int ihi = domainx.bigEnd(0);
511  int jlo = domainy.smallEnd(1);
512  int jhi = domainy.bigEnd(1);
513 
514 #ifdef _OPENMP
515 #pragma omp parallel if (Gpu::notInLaunchRegion())
516 #endif
517  for (MFIter mfi(S_cur_data[IntVars::cons],TilingIfNotGPU()); mfi.isValid(); ++mfi)
518  {
519  Box tbx = mfi.nodaltilebox(0);
520  Box tbx_lo, tbx_hi;
521  if (tbx.smallEnd(0) == ilo) {
522  tbx_lo = makeSlab(tbx,0,ilo);
523  }
524  if (tbx.bigEnd(0) == ihi) {
525  tbx_hi = makeSlab(tbx,0,ihi);
526  }
527 
528  Box tby = mfi.nodaltilebox(1);
529  Box tby_lo, tby_hi;
530  if (tby.smallEnd(1) == jlo) {
531  tby_lo = makeSlab(tby,1,jlo);
532  }
533  if (tby.bigEnd(1) == jhi) {
534  tby_hi = makeSlab(tby,1,jhi);
535  }
536 
537  Array4<Real> rhs_xmom = S_rhs[IntVars::xmom].array(mfi);
538  Array4<Real> rhs_ymom = S_rhs[IntVars::ymom].array(mfi);
539 
540  Array4<const Real> rhs_cons = S_rhs[IntVars::cons].const_array(mfi);
541  Array4<const Real> cons_arr = S_cur_data[IntVars::cons].const_array(mfi);
542 
543  const auto& bdatxlo_n = bdy_data_xlo[n_time ][ivarU].const_array();
544  const auto& bdatxlo_np1 = bdy_data_xlo[n_time_p1][ivarU].const_array();
545  const auto& bdatxhi_n = bdy_data_xhi[n_time ][ivarU].const_array();
546  const auto& bdatxhi_np1 = bdy_data_xhi[n_time_p1][ivarU].const_array();
547 
548  const auto& bdatylo_n = bdy_data_ylo[n_time ][ivarV].const_array();
549  const auto& bdatylo_np1 = bdy_data_ylo[n_time_p1][ivarV].const_array();
550  const auto& bdatyhi_n = bdy_data_yhi[n_time ][ivarV].const_array();
551  const auto& bdatyhi_np1 = bdy_data_yhi[n_time_p1][ivarV].const_array();
552 
553  ParallelFor(tbx_lo, tbx_hi,
554  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
555  {
556  Real rho_tend = rhs_cons(i,j,k);
557  Real rho_val = Real(0.5) * (cons_arr(i,j,k) + cons_arr(i-1,j,k));
558  Real u_tend, u_val;
559  if (btenxlo) {
560  u_tend = btenxlo(i,j,k,BCVars::xvel_bc);
561  u_val = bdatxlo(i,j,k,BCVars::xvel_bc);
562  } else {
563  u_tend = (bdatxlo_np1(i,j,k) - bdatxlo_n(i,j,k)) / bdy_time_interval;
564  u_val = oma * bdatxlo_n(i,j,k) + alpha * bdatxlo_np1(i,j,k);
565  }
566  rhs_xmom(i,j,k) = rho_val * u_tend + u_val * rho_tend;
567  },
568  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
569  {
570  Real rho_tend = rhs_cons(i,j,k);
571  Real rho_val = Real(0.5) * (cons_arr(i,j,k) + cons_arr(i-1,j,k));
572  Real u_tend, u_val;
573  if (btenxhi) {
574  u_tend = btenxhi(i,j,k,BCVars::xvel_bc);
575  u_val = bdatxhi(i,j,k,BCVars::xvel_bc);
576  } else {
577  u_tend = (bdatxhi_np1(i,j,k) - bdatxhi_n(i,j,k)) / bdy_time_interval;
578  u_val = oma * bdatxhi_n(i,j,k) + alpha * bdatxhi_np1(i,j,k);
579  }
580  rhs_xmom(i,j,k) = rho_val * u_tend + u_val * rho_tend;
581  });
582 
583  ParallelFor(tby_lo, tby_hi,
584  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
585  {
586  Real rho_tend = rhs_cons(i,j,k);
587  Real rho_val = Real(0.5) * (cons_arr(i,j,k) + cons_arr(i,j-1,k));
588  Real v_tend, v_val;
589  if (btenylo) {
590  v_tend = btenylo(i,j,k,BCVars::yvel_bc);
591  v_val = bdatylo(i,j,k,BCVars::yvel_bc);
592  } else {
593  v_tend = (bdatylo_np1(i,j,k) - bdatylo_n(i,j,k)) / bdy_time_interval;
594  v_val = oma * bdatylo_n(i,j,k) + alpha * bdatylo_np1(i,j,k);
595  }
596  rhs_ymom(i,j,k) = rho_val * v_tend + v_val * rho_tend;
597  },
598  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
599  {
600  Real rho_tend = rhs_cons(i,j,k);
601  Real rho_val = Real(0.5) * (cons_arr(i,j,k) + cons_arr(i,j-1,k));
602  Real v_tend, v_val;
603  if (btenyhi) {
604  v_tend = btenyhi(i,j,k,BCVars::yvel_bc);
605  v_val = bdatyhi(i,j,k,BCVars::yvel_bc);
606  } else {
607  v_tend = (bdatyhi_np1(i,j,k) - bdatyhi_n(i,j,k)) / bdy_time_interval;
608  v_val = oma * bdatyhi_n(i,j,k) + alpha * bdatyhi_np1(i,j,k);
609  }
610  rhs_ymom(i,j,k) = rho_val * v_tend + v_val * rho_tend;
611  });
612  } // mfi
613 }
RhoTheta_comp
#define RhoTheta_comp
Definition: ERF_IndexDefines.H:37
dx
const Real dx
Definition: ERF_InitCustomPert_ABL.H:23
rdOcp
const Real rdOcp
Definition: ERF_InitCustomPert_Bomex.H:16
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
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 &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 > &data_arr, const amrex::Array4< amrex::Real > &rhs_arr, const amrex::Real &c_p, const amrex::Real &rdOcp, const int bdy_moist_nudge_type=0)
Definition: ERF_Utils.H:236
BCVars::RhoTheta_bc_comp
@ RhoTheta_bc_comp
Definition: ERF_IndexDefines.H:88
BCVars::yvel_bc
@ yvel_bc
Definition: ERF_IndexDefines.H:103
BCVars::xvel_bc
@ xvel_bc
Definition: ERF_IndexDefines.H:102
RealBdyVars::U
@ U
Definition: ERF_IndexDefines.H:123
RealBdyVars::NumTypes
@ NumTypes
Definition: ERF_IndexDefines.H:127
RealBdyVars::T
@ T
Definition: ERF_IndexDefines.H:125
RealBdyVars::V
@ V
Definition: ERF_IndexDefines.H:124
Vars::xvel
@ xvel
Definition: ERF_IndexDefines.H:175
Vars::cons
@ cons
Definition: ERF_IndexDefines.H:174
Vars::yvel
@ yvel
Definition: ERF_IndexDefines.H:176
mp_wsm6::alpha
real(kind=kind_phys), parameter, public alpha
Definition: ERF_module_mp_wsm6.F90:44
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().