ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
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 int &set_width, 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 &bdy_time_interval, const Real &start_bdy_time, const Real &time, const Real &delta_t, int width, int set_width, const Geometry &geom, Vector< MultiFab > &S_rhs, Vector< MultiFab > &S_old_data, 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)
 
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 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
611 {
612  BL_PROFILE_REGION("fine_compute_interior_ghost_RHS()");
613 
614  // Relaxation constants
615  Real F1 = 1./(10.*delta_t);
616  Real F2 = 1./(50.*delta_t);
617 
618  // Vector of MFs to hold data (dm differs w/ fine patch)
619  Vector<MultiFab> fmf_p_v;
620 
621  // Loop over the variables
622  for (int ivar_idx = 0; ivar_idx < IntVars::NumTypes; ++ivar_idx)
623  {
624  // Fine mfs
625  MultiFab& fmf = S_data_f[ivar_idx];
626  MultiFab& rhs = S_rhs_f [ivar_idx];
627 
628  // NOTE: These temporary MFs and copy operations are horrible
629  // for memory usage and efficiency. However, we need to
630  // have access to ghost cells in the cons array to convert
631  // from primitive u/v/w to momentum. Furthermore, the BA
632  // for the fine patches in ERFFillPatcher don't match the
633  // BA for the data/RHS. For this reason, the data is copied
634  // to a vector of MFs (with ghost cells) so the BAs match
635  // the BA of data/RHS and we have access to rho to convert
636  // prim to conserved.
637 
638  // Temp MF on box (distribution map differs w/ fine patch)
639  int num_var = fmf.nComp();
640  fmf_p_v.emplace_back(fmf.boxArray(), fmf.DistributionMap(), num_var, fmf.nGrowVect());
641  MultiFab& fmf_p = fmf_p_v[ivar_idx];
642  MultiFab::Copy(fmf_p,fmf, 0, 0, num_var, fmf.nGrowVect());
643 
644  // Integer mask MF
645  int set_mask_val;
646  int relax_mask_val;
647  iMultiFab* mask;
648 
649  // Fill fine patch on interior halo region
650  //==========================================================
651  if (ivar_idx == IntVars::cons)
652  {
653  FPr_c->FillRelax(fmf_p, time, void_bc, domain_bcs_type);
654  mask = FPr_c->GetMask();
655  set_mask_val = FPr_c->GetSetMaskVal();
656  relax_mask_val = FPr_c->GetRelaxMaskVal();
657  }
658  else if (ivar_idx == IntVars::xmom)
659  {
660  FPr_u->FillRelax(fmf_p, time, void_bc, domain_bcs_type);
661  mask = FPr_u->GetMask();
662  set_mask_val = FPr_u->GetSetMaskVal();
663  relax_mask_val = FPr_u->GetRelaxMaskVal();
664 
665 #ifdef _OPENMP
666 #pragma omp parallel if (Gpu::notInLaunchRegion())
667 #endif
668  for ( MFIter mfi(fmf_p,TilingIfNotGPU()); mfi.isValid(); ++mfi)
669  {
670  Box vbx = mfi.validbox();
671 
672  const Array4<Real>& prim_arr = fmf_p.array(mfi);
673  const Array4<const Real>& rho_arr = fmf_p_v[0].const_array(mfi);
674  const Array4<const int>& mask_arr = mask->const_array(mfi);
675 
676  ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
677  {
678  if (mask_arr(i,j,k) == relax_mask_val) {
679  Real rho_interp = 0.5 * ( rho_arr(i-1,j,k) + rho_arr(i,j,k) );
680  prim_arr(i,j,k) *= rho_interp;
681  }
682  });
683  } // mfi
684  }
685  else if (ivar_idx == IntVars::ymom)
686  {
687  FPr_v->FillRelax(fmf_p, time, void_bc, domain_bcs_type);
688  mask = FPr_v->GetMask();
689  set_mask_val = FPr_v->GetSetMaskVal();
690  relax_mask_val = FPr_v->GetRelaxMaskVal();
691 
692 #ifdef _OPENMP
693 #pragma omp parallel if (Gpu::notInLaunchRegion())
694 #endif
695  for ( MFIter mfi(fmf_p,TilingIfNotGPU()); mfi.isValid(); ++mfi)
696  {
697  Box vbx = mfi.validbox();
698 
699  const Array4<Real>& prim_arr = fmf_p.array(mfi);
700  const Array4<const Real>& rho_arr = fmf_p_v[0].const_array(mfi);
701  const Array4<const int>& mask_arr = mask->const_array(mfi);
702 
703  ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
704  {
705  if (mask_arr(i,j,k) == relax_mask_val) {
706  Real rho_interp = 0.5 * ( rho_arr(i,j-1,k) + rho_arr(i,j,k) );
707  prim_arr(i,j,k) *= rho_interp;
708  }
709  });
710  } // mfi
711  }
712  else if (ivar_idx == IntVars::zmom)
713  {
714  FPr_w->FillRelax(fmf_p, time, void_bc, domain_bcs_type);
715  mask = FPr_w->GetMask();
716  set_mask_val = FPr_w->GetSetMaskVal();
717  relax_mask_val = FPr_w->GetRelaxMaskVal();
718 
719 #ifdef _OPENMP
720 #pragma omp parallel if (Gpu::notInLaunchRegion())
721 #endif
722  for ( MFIter mfi(fmf_p,TilingIfNotGPU()); mfi.isValid(); ++mfi)
723  {
724  Box vbx = mfi.validbox();
725 
726  const Array4<Real>& prim_arr = fmf_p.array(mfi);
727  const Array4<const Real>& rho_arr = fmf_p_v[0].const_array(mfi);
728  const Array4<const int>& mask_arr = mask->const_array(mfi);
729 
730  ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
731  {
732  if (mask_arr(i,j,k) == relax_mask_val) {
733  Real rho_interp = 0.5 * ( rho_arr(i,j,k-1) + rho_arr(i,j,k) );
734  prim_arr(i,j,k) *= rho_interp;
735  }
736  });
737  } // mfi
738  } else {
739  Abort("Dont recognize this variable type in fine_compute_interior_ghost_RHS");
740  }
741 
742 
743  // Zero RHS in set region
744  //==========================================================
745 #ifdef _OPENMP
746 #pragma omp parallel if (Gpu::notInLaunchRegion())
747 #endif
748  for ( MFIter mfi(rhs,TilingIfNotGPU()); mfi.isValid(); ++mfi)
749  {
750  Box vbx = mfi.validbox();
751  const Array4<Real>& rhs_arr = rhs.array(mfi);
752  const Array4<const int>& mask_arr = mask->const_array(mfi);
753 
754  ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
755  {
756  if (mask_arr(i,j,k) == set_mask_val) {
757  rhs_arr(i,j,k) = 0.0;
758  }
759  });
760  } // mfi
761 
762  // For Laplacian stencil
763  rhs.FillBoundary(geom.periodicity());
764 
765 
766  // Compute RHS in relaxation region
767  //==========================================================
768 #ifdef _OPENMP
769 #pragma omp parallel if (Gpu::notInLaunchRegion())
770 #endif
771  for ( MFIter mfi(fmf_p,TilingIfNotGPU()); mfi.isValid(); ++mfi)
772  {
773  Box vbx = mfi.validbox();
774  const Array4<Real>& rhs_arr = rhs.array(mfi);
775  const Array4<const Real>& fine_arr = fmf_p.const_array(mfi);
776  const Array4<const Real>& data_arr = fmf.const_array(mfi);
777  const Array4<const int>& mask_arr = mask->const_array(mfi);
778 
779  const auto& vbx_lo = lbound(vbx);
780  const auto& vbx_hi = ubound(vbx);
781 
782  int icomp = 0;
783 
784  int Spec_z = set_width;
785  int Relax_z = width - Spec_z;
786  Real num = Real(Spec_z + Relax_z);
787  Real denom = Real(Relax_z - 1);
788  ParallelFor(vbx, num_var, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
789  {
790  if (mask_arr(i,j,k) == relax_mask_val) {
791 
792  // Indices
793  Real n_ind(-1); // Set to -1 to quiet compiler warning
794  int ii(width-1); int jj(width-1);
795  bool near_x_lo_wall(false); bool near_x_hi_wall(false);
796  bool near_y_lo_wall(false); bool near_y_hi_wall(false);
797  bool mask_x_found(false); bool mask_y_found(false);
798 
799  // Near x-wall
800  if ((i-vbx_lo.x) < width) {
801  near_x_lo_wall = true;
802  ii = i-vbx_lo.x;
803  if (mask_arr(vbx_lo.x,j,k) == 2) mask_x_found = true;
804  } else if ((vbx_hi.x-i) < width) {
805  near_x_hi_wall = true;
806  ii = vbx_hi.x-i;
807  if (mask_arr(vbx_hi.x,j,k) == 2) mask_x_found = true;
808  }
809 
810  // Near y-wall
811  if ((j-vbx_lo.y) < width) {
812  near_y_lo_wall = true;
813  jj = j-vbx_lo.y;
814  if (mask_arr(i,vbx_lo.y,k) == 2) mask_y_found = true;
815  } else if ((vbx_hi.y-j) < width) {
816  near_y_hi_wall = true;
817  jj = vbx_hi.y-j;
818  if (mask_arr(i,vbx_hi.y,k) == 2) mask_y_found = true;
819  }
820 
821  // Found a nearby masked cell (valid n_ind)
822  if (mask_x_found && mask_y_found) {
823  n_ind = std::min(ii,jj) + 1.0;
824  } else if (mask_x_found) {
825  n_ind = ii + 1.0;
826  } else if (mask_y_found) {
827  n_ind = jj + 1.0;
828  // Pesky corner cell
829  } else {
830  if (near_x_lo_wall || near_x_hi_wall) {
831  Real dj_min{width-1.0};
832  int j_lb = std::max(vbx_lo.y,j-width);
833  int j_ub = std::min(vbx_hi.y,j+width);
834  int li = (near_x_lo_wall) ? vbx_lo.x : vbx_hi.x;
835  for (int lj(j_lb); lj<=j_ub; ++lj) {
836  if (mask_arr(li,lj,k) == 2) {
837  mask_y_found = true;
838  dj_min = std::min(dj_min,(Real) std::abs(lj-j));
839  }
840  }
841  if (mask_y_found) {
842  Real mag = sqrt( Real(dj_min*dj_min + ii*ii) );
843  n_ind = std::min(mag,width-1.0) + 1.0;
844  } else {
845  Abort("Mask not found near x wall!");
846  }
847  } else if (near_y_lo_wall || near_y_hi_wall) {
848  Real di_min{width-1.0};
849  int i_lb = std::max(vbx_lo.x,i-width);
850  int i_ub = std::min(vbx_hi.x,i+width);
851  int lj = (near_y_lo_wall) ? vbx_lo.y : vbx_hi.y;
852  for (int li(i_lb); li<=i_ub; ++li) {
853  if (mask_arr(li,lj,k) == 2) {
854  mask_x_found = true;
855  di_min = std::min(di_min,(Real) std::abs(li-i));
856  }
857  }
858  if (mask_x_found) {
859  Real mag = sqrt( Real(di_min*di_min + jj*jj) );
860  n_ind = std::min(mag,width-1.0) + 1.0;
861  } else {
862  Abort("Mask not found near y wall!");
863  }
864  } else {
865  Abort("Relaxation cell must be near a wall!");
866  }
867  }
868 
869  Real Factor = (num - n_ind)/denom;
870  Real d = data_arr(i ,j ,k ,n+icomp) + delta_t*rhs_arr(i , j , k ,n+icomp);
871  Real d_ip1 = data_arr(i+1,j ,k ,n+icomp) + delta_t*rhs_arr(i+1, j , k ,n+icomp);
872  Real d_im1 = data_arr(i-1,j ,k ,n+icomp) + delta_t*rhs_arr(i-1, j , k ,n+icomp);
873  Real d_jp1 = data_arr(i ,j+1,k ,n+icomp) + delta_t*rhs_arr(i , j+1, k ,n+icomp);
874  Real d_jm1 = data_arr(i ,j-1,k ,n+icomp) + delta_t*rhs_arr(i , j-1, k ,n+icomp);
875  Real delta = fine_arr(i ,j ,k,n) - d;
876  Real delta_xp = fine_arr(i+1,j ,k,n) - d_ip1;
877  Real delta_xm = fine_arr(i-1,j ,k,n) - d_im1;
878  Real delta_yp = fine_arr(i ,j+1,k,n) - d_jp1;
879  Real delta_ym = fine_arr(i ,j-1,k,n) - d_jm1;
880  Real Laplacian = delta_xp + delta_xm + delta_yp + delta_ym - 4.0*delta;
881  rhs_arr(i,j,k,n) += (F1*delta - F2*Laplacian) * Factor;
882  }
883  });
884  } // mfi
885  } // ivar_idx
886 }
PhysBCFunctNoOp void_bc
Definition: ERF_InteriorGhostCells.cpp:5
amrex::iMultiFab * GetMask()
Definition: ERF_FillPatcher.H:43
void FillRelax(amrex::MultiFab &mf, amrex::Real time, BC &cbc, amrex::Vector< amrex::BCRec > const &bcs)
Definition: ERF_FillPatcher.H:105
int GetSetMaskVal()
Definition: ERF_FillPatcher.H:39
int GetRelaxMaskVal()
Definition: ERF_FillPatcher.H:41
@ NumTypes
Definition: ERF_IndexDefines.H:154
@ ymom
Definition: ERF_IndexDefines.H:152
@ cons
Definition: ERF_IndexDefines.H:150
@ zmom
Definition: ERF_IndexDefines.H:153
@ xmom
Definition: ERF_IndexDefines.H:151
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
117 {
118  AMREX_ALWAYS_ASSERT(bx.ixType() == domain.ixType());
119 
120  // Domain bounds without ghost cells
121  const auto& dom_lo = lbound(domain);
122  const auto& dom_hi = ubound(domain);
123 
124  // Four boxes matching the domain
125  Box gdom_xlo(domain); Box gdom_xhi(domain);
126  Box gdom_ylo(domain); Box gdom_yhi(domain);
127 
128  // Get offsets from box index type
129  IntVect iv_type = bx.ixType().toIntVect();
130  int offx = (iv_type[0]==1) ? 0 : -1;
131  int offy = (iv_type[1]==1) ? 0 : -1;
132 
133  // Stagger the boxes based upon index type
134  gdom_xlo += IntVect(offx,0,0); gdom_xhi += IntVect(-offx,0,0);
135  gdom_ylo += IntVect(0,offy,0); gdom_yhi += IntVect(0,-offy,0);
136 
137  // Trim the boxes to only include internal ghost cells
138  gdom_xlo.setBig(0,dom_lo.x+set_width+offx-1); gdom_xhi.setSmall(0,dom_hi.x-set_width-offx+1);
139  gdom_ylo.setBig(1,dom_lo.y+set_width+offy-1); gdom_yhi.setSmall(1,dom_hi.y-set_width-offy+1);
140 
141  // Remove overlapping corners from y-face boxes
142  gdom_ylo.setSmall(0,gdom_xlo.bigEnd(0)+1); gdom_ylo.setBig(0,gdom_xhi.smallEnd(0)-1);
143  gdom_yhi.setSmall(0,gdom_xlo.bigEnd(0)+1); gdom_yhi.setBig(0,gdom_xhi.smallEnd(0)-1);
144 
145  // Grow boxes to get external ghost cells only
146  gdom_xlo.growLo(0,ng_vect[0]+offx); gdom_xhi.growHi(0,ng_vect[0]+offx);
147  gdom_xlo.grow (1,ng_vect[1] ); gdom_xhi.grow (1,ng_vect[1] );
148  gdom_ylo.growLo(1,ng_vect[1]+offy); gdom_yhi.growHi(1,ng_vect[1]+offy);
149 
150  // Populate everything
151  bx_xlo = (bx & gdom_xlo);
152  bx_xhi = (bx & gdom_xhi);
153  bx_ylo = (bx & gdom_ylo);
154  bx_yhi = (bx & gdom_yhi);
155 }

Referenced by realbdy_interior_bxs_xy().

Here is the caller graph for this function:

◆ realbdy_compute_interior_ghost_rhs()

void realbdy_compute_interior_ghost_rhs ( const Real &  bdy_time_interval,
const Real &  start_bdy_time,
const Real &  time,
const Real &  delta_t,
int  width,
int  set_width,
const Geometry &  geom,
Vector< MultiFab > &  S_rhs,
Vector< MultiFab > &  S_old_data,
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 
)

Compute the RHS in the relaxation zone

Parameters
[in]bdy_time_intervaltime interval between boundary condition time stamps
[in]timecurrent time
[in]delta_ttimestep
[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
[in]start_bdy_timetime of the first boundary data read in
189 {
190  BL_PROFILE_REGION("realbdy_compute_interior_ghost_RHS()");
191 
192  // NOTE: We pass the full width into this routine.
193  // For relaxation, the last cell is a halo
194  // cell for the Laplacian. We remove that
195  // cell here if it is present.
196 
197  // The width to do RHS augmentation
198  if (width > set_width+1) width -= 1;
199 
200  // Relaxation constants
201  Real F1 = 1./(10.*delta_t);
202  Real F2 = 1./(50.*delta_t);
203 
204  // Time interpolation
205  Real dT = bdy_time_interval;
206  Real time_since_start = time - start_bdy_time;
207  int n_time = static_cast<int>( time_since_start / dT);
208  Real alpha = (time_since_start - n_time * dT) / dT;
209  AMREX_ALWAYS_ASSERT( alpha >= 0. && alpha <= 1.0);
210  Real oma = 1.0 - alpha;
211 
212  /*
213  // UNIT TEST DEBUG
214  oma = 1.0; alpha = 0.0;
215  */
216 
217  // Temporary FABs for storage (owned/filled on all ranks)
218  FArrayBox U_xlo, U_xhi, U_ylo, U_yhi;
219  FArrayBox V_xlo, V_xhi, V_ylo, V_yhi;
220  FArrayBox T_xlo, T_xhi, T_ylo, T_yhi;
221 
222  // Variable index map (WRFBdyVars -> Vars)
223  Vector<int> var_map = {Vars::xvel, Vars::yvel, Vars::cons, Vars::cons};
224  Vector<int> ivar_map = {IntVars::xmom, IntVars::ymom, IntVars::cons, IntVars::cons};
225 
226  // Variable icomp map
227  Vector<int> comp_map = {0, 0, RhoTheta_comp};
228 
229  // Indices
230  int ivarU = RealBdyVars::U;
231  int ivarV = RealBdyVars::V;
232  int ivarT = RealBdyVars::T;
233  int BdyEnd = RealBdyVars::NumTypes-1;
234 
235 
236  // NOTE: These sizing of the temporary BDY FABS is
237  // GLOBAL and occurs over the entire BDY region.
238 
239  // Size the FABs
240  //==========================================================
241  for (int ivar(ivarU); ivar < BdyEnd; ivar++) {
242  int ivar_idx = var_map[ivar];
243  Box domain = geom.Domain();
244  auto ixtype = S_cur_data[ivar_idx].boxArray().ixType();
245  domain.convert(ixtype);
246 
247  // Grown domain to get the 4 halo boxes w/ ghost cells
248  // NOTE: 2 ghost cells needed here for Laplacian
249  // halo cell.
250  IntVect ng_vect{2,2,0};
251  Box gdom(domain); gdom.grow(ng_vect);
252  Box bx_xlo, bx_xhi, bx_ylo, bx_yhi;
253  realbdy_interior_bxs_xy(gdom, domain, width,
254  bx_xlo, bx_xhi,
255  bx_ylo, bx_yhi,
256  0, ng_vect, true);
257 
258  // Size the FABs
259  if (ivar == ivarU) {
260  U_xlo.resize(bx_xlo,1,The_Async_Arena()); U_xhi.resize(bx_xhi,1,The_Async_Arena());
261  U_ylo.resize(bx_ylo,1,The_Async_Arena()); U_yhi.resize(bx_yhi,1,The_Async_Arena());
262  } else if (ivar == ivarV) {
263  V_xlo.resize(bx_xlo,1,The_Async_Arena()); V_xhi.resize(bx_xhi,1,The_Async_Arena());
264  V_ylo.resize(bx_ylo,1,The_Async_Arena()); V_yhi.resize(bx_yhi,1,The_Async_Arena());
265  } else if (ivar == ivarT){
266  T_xlo.resize(bx_xlo,1,The_Async_Arena()); T_xhi.resize(bx_xhi,1,The_Async_Arena());
267  T_ylo.resize(bx_ylo,1,The_Async_Arena()); T_yhi.resize(bx_yhi,1,The_Async_Arena());
268  } else {
269  continue;
270  }
271  } // ivar
272 
273 
274  // NOTE: These operations use the BDY FABS and RHO. The
275  // use of RHO to go from PRIM -> CONS requires that
276  // these operations be LOCAL. So we have allocated
277  // enough space to do global operations (1 rank) but
278  // will fill a subset of that data that the rank owns.
279 
280  // Populate FABs from bdy interpolation (primitive vars)
281  //==========================================================
282  for (int ivar(ivarU); ivar < BdyEnd; ivar++) {
283  int ivar_idx = var_map[ivar];
284  Box domain = geom.Domain();
285  auto ixtype = S_cur_data[ivar_idx].boxArray().ixType();
286  domain.convert(ixtype);
287  const auto& dom_lo = lbound(domain);
288  const auto& dom_hi = ubound(domain);
289 
290 #ifdef _OPENMP
291 #pragma omp parallel if (Gpu::notInLaunchRegion())
292 #endif
293  for (MFIter mfi(S_cur_data[ivar_idx],TilingIfNotGPU()); mfi.isValid(); ++mfi) {
294  // We need lateral ghost cells for the Laplacian
295  // NOTE: We don't write into the ghost cells
296  IntVect ng_vect{2,2,0};
297  Box gtbx = grow(mfi.tilebox(ixtype.toIntVect()),ng_vect);
298  Box tbx_xlo, tbx_xhi, tbx_ylo, tbx_yhi;
299  realbdy_interior_bxs_xy(gtbx, domain, width,
300  tbx_xlo, tbx_xhi,
301  tbx_ylo, tbx_yhi,
302  0, ng_vect, true);
303 
304  Array4<Real> arr_xlo; Array4<Real> arr_xhi;
305  Array4<Real> arr_ylo; Array4<Real> arr_yhi;
306  if (ivar == ivarU) {
307  arr_xlo = U_xlo.array(); arr_xhi = U_xhi.array();
308  arr_ylo = U_ylo.array(); arr_yhi = U_yhi.array();
309  } else if (ivar == ivarV) {
310  arr_xlo = V_xlo.array(); arr_xhi = V_xhi.array();
311  arr_ylo = V_ylo.array(); arr_yhi = V_yhi.array();
312  } else if (ivar == ivarT){
313  arr_xlo = T_xlo.array(); arr_xhi = T_xhi.array();
314  arr_ylo = T_ylo.array(); arr_yhi = T_yhi.array();
315  } else {
316  continue;
317  }
318 
319  // Boundary data at fixed time intervals
320  const auto& bdatxlo_n = bdy_data_xlo[n_time ][ivar].const_array();
321  const auto& bdatxlo_np1 = bdy_data_xlo[n_time+1][ivar].const_array();
322  const auto& bdatxhi_n = bdy_data_xhi[n_time ][ivar].const_array();
323  const auto& bdatxhi_np1 = bdy_data_xhi[n_time+1][ivar].const_array();
324  const auto& bdatylo_n = bdy_data_ylo[n_time ][ivar].const_array();
325  const auto& bdatylo_np1 = bdy_data_ylo[n_time+1][ivar].const_array();
326  const auto& bdatyhi_n = bdy_data_yhi[n_time ][ivar].const_array();
327  const auto& bdatyhi_np1 = bdy_data_yhi[n_time+1][ivar].const_array();
328 
329  // Current density to convert to conserved vars
330  Array4<Real> r_arr = S_cur_data[IntVars::cons].array(mfi);
331 
332  // NOTE: width is now one less than the total bndy width
333  // if we have a relaxation zone; so we can access
334  // dom_lo/hi +- width. If we do not have a relax
335  // zone, this offset is set_width - 1.
336  int offset = set_width - 1;
337  if (width > set_width) offset = width;
338 
339  // Populate with interpolation (protect from ghost cells)
340  ParallelFor(tbx_xlo, tbx_xhi,
341  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
342  {
343  int ii = std::max(i , dom_lo.x);
344  ii = std::min(ii, dom_lo.x+offset);
345  int jj = std::max(j , dom_lo.y);
346  jj = std::min(jj, dom_hi.y);
347 
348  Real rho_interp;
349  if (ivar==ivarU) {
350  rho_interp = 0.5 * ( r_arr(i-1,j ,k) + r_arr(i,j,k) );
351  } else if (ivar==ivarV) {
352  rho_interp = 0.5 * ( r_arr(i ,j-1,k) + r_arr(i,j,k) );
353  } else {
354  rho_interp = r_arr(i,j,k);
355  }
356  arr_xlo(i,j,k) = rho_interp * ( oma * bdatxlo_n (ii,jj,k,0)
357  + alpha * bdatxlo_np1(ii,jj,k,0) );
358  },
359  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
360  {
361  int ii = std::max(i , dom_hi.x-offset);
362  ii = std::min(ii, dom_hi.x);
363  int jj = std::max(j , dom_lo.y);
364  jj = std::min(jj, dom_hi.y);
365 
366  Real rho_interp;
367  if (ivar==ivarU) {
368  rho_interp = 0.5 * ( r_arr(i-1,j ,k) + r_arr(i,j,k) );
369  } else if (ivar==ivarV) {
370  rho_interp = 0.5 * ( r_arr(i ,j-1,k) + r_arr(i,j,k) );
371  } else {
372  rho_interp = r_arr(i,j,k);
373  }
374  arr_xhi(i,j,k) = rho_interp * ( oma * bdatxhi_n (ii,jj,k,0)
375  + alpha * bdatxhi_np1(ii,jj,k,0) );
376  });
377 
378  ParallelFor(tbx_ylo, tbx_yhi,
379  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
380  {
381  int ii = std::max(i , dom_lo.x);
382  ii = std::min(ii, dom_hi.x);
383  int jj = std::max(j , dom_lo.y);
384  jj = std::min(jj, dom_lo.y+offset);
385 
386  Real rho_interp;
387  if (ivar==ivarU) {
388  rho_interp = 0.5 * ( r_arr(i-1,j ,k) + r_arr(i,j,k) );
389  } else if (ivar==ivarV) {
390  rho_interp = 0.5 * ( r_arr(i ,j-1,k) + r_arr(i,j,k) );
391  } else {
392  rho_interp = r_arr(i,j,k);
393  }
394  arr_ylo(i,j,k) = rho_interp * ( oma * bdatylo_n (ii,jj,k,0)
395  + alpha * bdatylo_np1(ii,jj,k,0) );
396  },
397  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
398  {
399  int ii = std::max(i , dom_lo.x);
400  ii = std::min(ii, dom_hi.x);
401  int jj = std::max(j , dom_hi.y-offset);
402  jj = std::min(jj, dom_hi.y);
403 
404  Real rho_interp;
405  if (ivar==ivarU) {
406  rho_interp = 0.5 * ( r_arr(i-1,j ,k) + r_arr(i,j,k) );
407  } else if (ivar==ivarV) {
408  rho_interp = 0.5 * ( r_arr(i ,j-1,k) + r_arr(i,j,k) );
409  } else {
410  rho_interp = r_arr(i,j,k);
411  }
412  arr_yhi(i,j,k) = rho_interp * ( oma * bdatyhi_n (ii,jj,k,0)
413  + alpha * bdatyhi_np1(ii,jj,k,0) );
414  });
415  } // mfi
416  } // ivar
417 
418 
419  // NOTE: These operations use current RHS, so they are
420  // LOCAL and occur over the data owned by a given rank.
421 
422  // Compute RHS in specified region
423  //==========================================================
424  if (set_width > 0) {
425  for (int ivar(ivarU); ivar < BdyEnd; ivar++) {
426  int ivar_idx = ivar_map[ivar];
427  int icomp = comp_map[ivar];
428 
429  Box domain = geom.Domain();
430  auto ix_type = S_old_data[ivar_idx].boxArray().ixType();
431  auto iv_type = ix_type.toIntVect();
432  domain.convert(ix_type);
433  const auto& dom_hi = ubound(domain);
434  const auto& dom_lo = lbound(domain);
435 
436  int set_width_x = (iv_type[0]) ? set_width : set_width-1;
437  int set_width_y = (iv_type[1]) ? set_width : set_width-1;
438 
439 #ifdef _OPENMP
440 #pragma omp parallel if (Gpu::notInLaunchRegion())
441 #endif
442  for (MFIter mfi(S_old_data[ivar_idx],TilingIfNotGPU()); mfi.isValid(); ++mfi) {
443  Box tbx = mfi.tilebox();
444  Box tbx_xlo, tbx_xhi, tbx_ylo, tbx_yhi;
445  realbdy_interior_bxs_xy(tbx, domain, width,
446  tbx_xlo, tbx_xhi,
447  tbx_ylo, tbx_yhi);
448 
449  Array4<Real> rhs_arr; Array4<Real> data_arr;
450  Array4<Real> arr_xlo; Array4<Real> arr_xhi;
451  Array4<Real> arr_ylo; Array4<Real> arr_yhi;
452  if (ivar == ivarU) {
453  arr_xlo = U_xlo.array(); arr_xhi = U_xhi.array();
454  arr_ylo = U_ylo.array(); arr_yhi = U_yhi.array();
455  rhs_arr = S_rhs[IntVars::xmom].array(mfi);
456  data_arr = S_old_data[IntVars::xmom].array(mfi);
457  } else if (ivar == ivarV) {
458  arr_xlo = V_xlo.array(); arr_xhi = V_xhi.array();
459  arr_ylo = V_ylo.array(); arr_yhi = V_yhi.array();
460  rhs_arr = S_rhs[IntVars::ymom].array(mfi);
461  data_arr = S_old_data[IntVars::ymom].array(mfi);
462  } else if (ivar == ivarT){
463  arr_xlo = T_xlo.array(); arr_xhi = T_xhi.array();
464  arr_ylo = T_ylo.array(); arr_yhi = T_yhi.array();
465  rhs_arr = S_rhs[IntVars::cons].array(mfi);
466  data_arr = S_old_data[IntVars::cons].array(mfi);
467  } else {
468  continue;
469  }
470 
471  realbdy_set_rhs_in_spec_region(delta_t, icomp, 1,
472  width, set_width_x, set_width_y,
473  dom_lo, dom_hi,
474  tbx_xlo, tbx_xhi, tbx_ylo, tbx_yhi,
475  arr_xlo, arr_xhi, arr_ylo, arr_yhi,
476  data_arr, rhs_arr);
477 
478  } // mfi
479  } // ivar
480  } // set_width
481 
482  // NOTE: These operations use current density, so they are
483  // LOCAL and occur over the data owned by a given rank
484 
485  // Compute RHS in relaxation region
486  //==========================================================
487  if (width > set_width) {
488  for (int ivar(ivarU); ivar < BdyEnd; ivar++) {
489  int ivar_idx = ivar_map[ivar];
490  int icomp = comp_map[ivar];
491 
492  Box domain = geom.Domain();
493  domain.convert(S_cur_data[ivar_idx].boxArray().ixType());
494  const auto& dom_hi = ubound(domain);
495  const auto& dom_lo = lbound(domain);
496  IntVect ng_vect = S_cur_data[ivar_idx].nGrowVect();
497 
498 #ifdef _OPENMP
499 #pragma omp parallel if (Gpu::notInLaunchRegion())
500 #endif
501  for (MFIter mfi(S_cur_data[ivar_idx],TilingIfNotGPU()); mfi.isValid(); ++mfi) {
502  Box tbx = mfi.tilebox();
503  Box tbx_xlo, tbx_xhi, tbx_ylo, tbx_yhi;
504  realbdy_interior_bxs_xy(tbx, domain, width,
505  tbx_xlo, tbx_xhi,
506  tbx_ylo, tbx_yhi,
507  set_width, ng_vect);
508 
509  Array4<Real> rhs_arr; Array4<Real> data_arr;
510  Array4<Real> arr_xlo; Array4<Real> arr_xhi;
511  Array4<Real> arr_ylo; Array4<Real> arr_yhi;
512  if (ivar == ivarU) {
513  arr_xlo = U_xlo.array(); arr_xhi = U_xhi.array();
514  arr_ylo = U_ylo.array(); arr_yhi = U_yhi.array();
515  rhs_arr = S_rhs[IntVars::xmom].array(mfi);
516  data_arr = S_cur_data[IntVars::xmom].array(mfi);
517  } else if (ivar == ivarV) {
518  arr_xlo = V_xlo.array(); arr_xhi = V_xhi.array();
519  arr_ylo = V_ylo.array(); arr_yhi = V_yhi.array();
520  rhs_arr = S_rhs[IntVars::ymom].array(mfi);
521  data_arr = S_cur_data[IntVars::ymom].array(mfi);
522  } else if (ivar == ivarT){
523  arr_xlo = T_xlo.array(); arr_xhi = T_xhi.array();
524  arr_ylo = T_ylo.array(); arr_yhi = T_yhi.array();
525  rhs_arr = S_rhs[IntVars::cons].array(mfi);
526  data_arr = S_cur_data[IntVars::cons].array(mfi);
527  } else {
528  continue;
529  }
530 
532  width, set_width, dom_lo, dom_hi, F1, F2,
533  tbx_xlo, tbx_xhi, tbx_ylo, tbx_yhi,
534  arr_xlo, arr_xhi, arr_ylo, arr_yhi,
535  data_arr, rhs_arr);
536 
537  /*
538  // UNIT TEST DEBUG
539  realbdy_interior_bxs_xy(tbx, domain, width+1,
540  tbx_xlo, tbx_xhi,
541  tbx_ylo, tbx_yhi);
542  ParallelFor(tbx_xlo, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
543  {
544  if (arr_xlo(i,j,k) != data_arr(i,j,k,icomp)) {
545  Print() << "ERROR XLO: " << ivar << ' ' << icomp << ' ' << IntVect(i,j,k) << "\n";
546  Print() << "DATA: " << data_arr(i,j,k,icomp) << ' ' << arr_xlo(i,j,k) << "\n";
547  exit(0);
548  }
549  });
550  ParallelFor(tbx_xhi, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
551  {
552  if (arr_xhi(i,j,k) != data_arr(i,j,k,icomp)) {
553  Print() << "ERROR XHI: " << ivar << ' ' << icomp << ' ' << IntVect(i,j,k) << "\n";
554  Print() << "DATA: " << data_arr(i,j,k,icomp) << ' ' << arr_xhi(i,j,k) << "\n";
555  exit(0);
556  }
557  });
558  ParallelFor(tbx_ylo, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
559  {
560  if (arr_ylo(i,j,k) != data_arr(i,j,k,icomp)) {
561  Print() << "ERROR YLO: " << ivar << ' ' << icomp << ' ' << IntVect(i,j,k) << "\n";
562  Print() << "DATA: " << data_arr(i,j,k,icomp) << ' ' << arr_ylo(i,j,k) << "\n";
563  exit(0);
564  }
565  });
566  ParallelFor(tbx_yhi, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
567  {
568  if (arr_yhi(i,j,k) != data_arr(i,j,k,icomp)) {
569  Print() << "ERROR YHI: " << ivar << ' ' << icomp << ' ' << IntVect(i,j,k) << "\n";
570  Print() << "DATA: " << data_arr(i,j,k,icomp) << ' ' << arr_yhi(i,j,k) << "\n";
571  exit(0);
572  }
573  });
574  */
575  } // mfi
576  } // ivar
577  } // width
578  //ParallelDescriptor::Barrier();
579  //exit(0);
580 }
#define RhoTheta_comp
Definition: ERF_IndexDefines.H:37
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 int &set_width, const IntVect &ng_vect, const bool get_int_ng)
Definition: ERF_InteriorGhostCells.cpp:23
AMREX_FORCE_INLINE IntVect offset(const int face_dir, const int normal)
Definition: ERF_ReadBndryPlanes.cpp:28
AMREX_GPU_HOST AMREX_FORCE_INLINE void realbdy_set_rhs_in_spec_region(const amrex::Real &dt, const int &icomp, const int &num_var, const int &width, const int &set_width_x, const int &set_width_y, const amrex::Dim3 &dom_lo, const amrex::Dim3 &dom_hi, 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)
Definition: ERF_Utils.H:169
AMREX_GPU_HOST AMREX_FORCE_INLINE void realbdy_compute_laplacian_relaxation(const int &icomp, const int &num_var, const int &width, const int &set_width, const amrex::Dim3 &dom_lo, const amrex::Dim3 &dom_hi, const amrex::Real &F1, const amrex::Real &F2, 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)
Definition: ERF_Utils.H:263
@ U
Definition: ERF_IndexDefines.H:108
@ NumTypes
Definition: ERF_IndexDefines.H:112
@ T
Definition: ERF_IndexDefines.H:110
@ V
Definition: ERF_IndexDefines.H:109
@ xvel
Definition: ERF_IndexDefines.H:141
@ cons
Definition: ERF_IndexDefines.H:140
@ yvel
Definition: ERF_IndexDefines.H:142
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 int &  set_width,
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
33 {
34  AMREX_ALWAYS_ASSERT(bx.ixType() == domain.ixType());
35 
36  //==================================================================
37  // NOTE: X-face boxes take ownership of the overlapping region.
38  // With exterior ghost cells (ng_vect != 0), the x-face
39  // boxes will have exterior ghost cells in both x & y.
40  //==================================================================
41 
42  // Domain bounds without ghost cells
43  const auto& dom_lo = lbound(domain);
44  const auto& dom_hi = ubound(domain);
45 
46  // Four boxes matching the domain
47  Box gdom_xlo(domain); Box gdom_xhi(domain);
48  Box gdom_ylo(domain); Box gdom_yhi(domain);
49 
50  // Trim the boxes to only include internal ghost cells
51  gdom_xlo.setBig(0,dom_lo.x+width-1); gdom_xhi.setSmall(0,dom_hi.x-width+1);
52  gdom_ylo.setBig(1,dom_lo.y+width-1); gdom_yhi.setSmall(1,dom_hi.y-width+1);
53 
54  // Remove overlapping corners from y-face boxes
55  gdom_ylo.setSmall(0,gdom_xlo.bigEnd(0)+1); gdom_ylo.setBig(0,gdom_xhi.smallEnd(0)-1);
56  gdom_yhi.setSmall(0,gdom_xlo.bigEnd(0)+1); gdom_yhi.setBig(0,gdom_xhi.smallEnd(0)-1);
57 
58  if (set_width>0) {
59  // Cut out the set region (uses offsets)
60  IntVect iv_type = bx.ixType().toIntVect();
61  Box sdom_xlo,sdom_xhi,sdom_ylo,sdom_yhi;
62  realbdy_bc_bxs_xy (grow(domain,ng_vect), domain, set_width,
63  sdom_xlo, sdom_xhi,
64  sdom_ylo, sdom_yhi,
65  ng_vect);
66  gdom_xlo.setSmall(0,sdom_xlo.bigEnd(0)+1); gdom_xhi.setBig(0,sdom_xhi.smallEnd(0)-1);
67  gdom_ylo.setSmall(1,sdom_ylo.bigEnd(1)+1); gdom_yhi.setBig(1,sdom_yhi.smallEnd(1)-1);
68  if (iv_type[1]==1) {
69  gdom_xlo.setSmall(1,sdom_ylo.bigEnd(1)+1); gdom_xhi.setSmall(1,sdom_ylo.bigEnd(1)+1);
70  gdom_xlo.setBig(1,sdom_yhi.smallEnd(1)-1); gdom_xhi.setBig(1,sdom_yhi.smallEnd(1)-1);
71  }
72  } else {
73  // Grow boxes to get external ghost cells only
74  gdom_xlo.growLo(0,ng_vect[0]); gdom_xhi.growHi(0,ng_vect[0]);
75  gdom_xlo.grow (1,ng_vect[1]); gdom_xhi.grow (1,ng_vect[1]);
76  gdom_ylo.growLo(1,ng_vect[1]); gdom_yhi.growHi(1,ng_vect[1]);
77  }
78 
79  // Grow boxes to get internal ghost cells
80  if (get_int_ng) {
81  gdom_xlo.growHi(0,ng_vect[0]); gdom_xhi.growLo(0,ng_vect[0]);
82  gdom_ylo.grow (0,ng_vect[0]); gdom_yhi.grow (0,ng_vect[0]);
83  gdom_ylo.growHi(1,ng_vect[1]); gdom_yhi.growLo(1,ng_vect[1]);
84  }
85 
86  // Populate everything
87  bx_xlo = (bx & gdom_xlo);
88  bx_xhi = (bx & gdom_xhi);
89  bx_ylo = (bx & gdom_ylo);
90  bx_yhi = (bx & gdom_yhi);
91 }
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)
Definition: ERF_InteriorGhostCells.cpp:109

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