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 double &time, const double &delta_t_d, const double &start_bdy_time, const double &final_bdy_time, const double &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, const bool use_wrf_bdy_density, const Real &bdy_rho_nudge_factor)
 
void fine_compute_interior_ghost_rhs (const double &time, const double &delta_t_d, 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 double &  time,
const double &  delta_t_d,
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
765 {
766  BL_PROFILE_REGION("fine_compute_interior_ghost_RHS()");
767 
768  Real delta_t = static_cast<Real>(delta_t_d);
769 
770  // Relaxation constants
771  Real F1 = one/(Real(10.)*delta_t);
772  Real F2 = one/(Real(50.)*delta_t);
773 
774  // Vector of MFs to hold data (dm differs w/ fine patch)
775  Vector<MultiFab> fmf_p_v;
776 
777  // Loop over the variables
778  for (int ivar_idx = 0; ivar_idx < IntVars::NumTypes; ++ivar_idx)
779  {
780  // Fine mfs
781  MultiFab& fmf = S_data_f[ivar_idx];
782  MultiFab& rhs = S_rhs_f [ivar_idx];
783 
784  // NOTE: These temporary MFs and copy operations are horrible
785  // for memory usage and efficiency. However, we need to
786  // have access to ghost cells in the cons array to convert
787  // from primitive u/v/w to momentum. Furthermore, the BA
788  // for the fine patches in ERFFillPatcher don't match the
789  // BA for the data/RHS. For this reason, the data is copied
790  // to a vector of MFs (with ghost cells) so the BAs match
791  // the BA of data/RHS and we have access to rho to convert
792  // prim to conserved.
793 
794  // Temp MF on box (distribution map differs w/ fine patch)
795  int num_var = fmf.nComp();
796  fmf_p_v.emplace_back(fmf.boxArray(), fmf.DistributionMap(), num_var, fmf.nGrowVect());
797  MultiFab& fmf_p = fmf_p_v[ivar_idx];
798  MultiFab::Copy(fmf_p,fmf, 0, 0, num_var, fmf.nGrowVect());
799 
800  // Integer mask MF
801  int set_mask_val;
802  int relax_mask_val;
803  iMultiFab* mask;
804 
805  // Fill fine patch on interior halo region
806  //==========================================================
807  if (ivar_idx == IntVars::cons)
808  {
809  FPr_c->FillRelax(fmf_p, time, void_bc, domain_bcs_type);
810  mask = FPr_c->GetMask();
811  set_mask_val = FPr_c->GetSetMaskVal();
812  relax_mask_val = FPr_c->GetRelaxMaskVal();
813  }
814  else if (ivar_idx == IntVars::xmom)
815  {
816  FPr_u->FillRelax(fmf_p, time, void_bc, domain_bcs_type);
817  mask = FPr_u->GetMask();
818  set_mask_val = FPr_u->GetSetMaskVal();
819  relax_mask_val = FPr_u->GetRelaxMaskVal();
820 
821 #ifdef _OPENMP
822 #pragma omp parallel if (Gpu::notInLaunchRegion())
823 #endif
824  for ( MFIter mfi(fmf_p,TilingIfNotGPU()); mfi.isValid(); ++mfi)
825  {
826  Box tbx = mfi.tilebox();
827 
828  const Array4<Real>& prim_arr = fmf_p.array(mfi);
829  const Array4<const Real>& rho_arr = fmf_p_v[0].const_array(mfi);
830  const Array4<const int>& mask_arr = mask->const_array(mfi);
831 
832  ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
833  {
834  if (mask_arr(i,j,k) == relax_mask_val) {
835  Real rho_interp = myhalf * ( rho_arr(i-1,j,k) + rho_arr(i,j,k) );
836  prim_arr(i,j,k) *= rho_interp;
837  }
838  });
839  } // mfi
840  }
841  else if (ivar_idx == IntVars::ymom)
842  {
843  FPr_v->FillRelax(fmf_p, time, void_bc, domain_bcs_type);
844  mask = FPr_v->GetMask();
845  set_mask_val = FPr_v->GetSetMaskVal();
846  relax_mask_val = FPr_v->GetRelaxMaskVal();
847 
848 #ifdef _OPENMP
849 #pragma omp parallel if (Gpu::notInLaunchRegion())
850 #endif
851  for ( MFIter mfi(fmf_p,TilingIfNotGPU()); mfi.isValid(); ++mfi)
852  {
853  Box tbx = mfi.tilebox();
854 
855  const Array4<Real>& prim_arr = fmf_p.array(mfi);
856  const Array4<const Real>& rho_arr = fmf_p_v[0].const_array(mfi);
857  const Array4<const int>& mask_arr = mask->const_array(mfi);
858 
859  ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
860  {
861  if (mask_arr(i,j,k) == relax_mask_val) {
862  Real rho_interp = myhalf * ( rho_arr(i,j-1,k) + rho_arr(i,j,k) );
863  prim_arr(i,j,k) *= rho_interp;
864  }
865  });
866  } // mfi
867  }
868  else if (ivar_idx == IntVars::zmom)
869  {
870  FPr_w->FillRelax(fmf_p, time, void_bc, domain_bcs_type);
871  mask = FPr_w->GetMask();
872  set_mask_val = FPr_w->GetSetMaskVal();
873  relax_mask_val = FPr_w->GetRelaxMaskVal();
874 
875 #ifdef _OPENMP
876 #pragma omp parallel if (Gpu::notInLaunchRegion())
877 #endif
878  for ( MFIter mfi(fmf_p,TilingIfNotGPU()); mfi.isValid(); ++mfi)
879  {
880  Box tbx = mfi.tilebox();
881 
882  const Array4<Real>& prim_arr = fmf_p.array(mfi);
883  const Array4<const Real>& rho_arr = fmf_p_v[0].const_array(mfi);
884  const Array4<const int>& mask_arr = mask->const_array(mfi);
885 
886  ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
887  {
888  if (mask_arr(i,j,k) == relax_mask_val) {
889  Real rho_interp = myhalf * ( rho_arr(i,j,k-1) + rho_arr(i,j,k) );
890  prim_arr(i,j,k) *= rho_interp;
891  }
892  });
893  } // mfi
894  } else {
895  Abort("Dont recognize this variable type in fine_compute_interior_ghost_RHS");
896  }
897 
898 
899  // Zero RHS in set region
900  //==========================================================
901 #ifdef _OPENMP
902 #pragma omp parallel if (Gpu::notInLaunchRegion())
903 #endif
904  for ( MFIter mfi(rhs,TilingIfNotGPU()); mfi.isValid(); ++mfi)
905  {
906  Box tbx = mfi.tilebox();
907  const Array4<Real>& rhs_arr = rhs.array(mfi);
908  const Array4<const int>& mask_arr = mask->const_array(mfi);
909 
910  ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
911  {
912  if (mask_arr(i,j,k) == set_mask_val) {
913  rhs_arr(i,j,k) = zero;
914  }
915  });
916  } // mfi
917 
918  // For Laplacian stencil
919  rhs.FillBoundary(geom.periodicity());
920 
921 
922  // Compute RHS in relaxation region
923  //==========================================================
924 #ifdef _OPENMP
925 #pragma omp parallel if (Gpu::notInLaunchRegion())
926 #endif
927  for ( MFIter mfi(fmf_p,TilingIfNotGPU()); mfi.isValid(); ++mfi)
928  {
929  Box tbx = mfi.tilebox();
930  const Array4<Real>& rhs_arr = rhs.array(mfi);
931  const Array4<const Real>& fine_arr = fmf_p.const_array(mfi);
932  const Array4<const Real>& data_arr = fmf.const_array(mfi);
933  const Array4<const int>& mask_arr = mask->const_array(mfi);
934 
935  Box vbx = mfi.validbox();
936  const auto& vbx_lo = lbound(vbx);
937  const auto& vbx_hi = ubound(vbx);
938 
939  int icomp = 0;
940 
941  int Spec_z = set_width;
942  int Relax_z = width - Spec_z;
943  Real num = Real(Spec_z + Relax_z);
944  Real denom = Real(Relax_z - 1);
945  ParallelFor(tbx, num_var, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
946  {
947  if (mask_arr(i,j,k) == relax_mask_val) {
948 
949  // Indices
950  Real n_ind(-1); // Set to -1 to quiet compiler warning
951  int ii(width-1); int jj(width-1);
952  bool near_x_lo_wall(false); bool near_x_hi_wall(false);
953  bool near_y_lo_wall(false); bool near_y_hi_wall(false);
954  bool mask_x_found(false); bool mask_y_found(false);
955 
956  // Near x-wall
957  if ((i-vbx_lo.x) < width) {
958  near_x_lo_wall = true;
959  ii = i-vbx_lo.x;
960  if (mask_arr(vbx_lo.x,j,k) == 2) mask_x_found = true;
961  } else if ((vbx_hi.x-i) < width) {
962  near_x_hi_wall = true;
963  ii = vbx_hi.x-i;
964  if (mask_arr(vbx_hi.x,j,k) == 2) mask_x_found = true;
965  }
966 
967  // Near y-wall
968  if ((j-vbx_lo.y) < width) {
969  near_y_lo_wall = true;
970  jj = j-vbx_lo.y;
971  if (mask_arr(i,vbx_lo.y,k) == 2) mask_y_found = true;
972  } else if ((vbx_hi.y-j) < width) {
973  near_y_hi_wall = true;
974  jj = vbx_hi.y-j;
975  if (mask_arr(i,vbx_hi.y,k) == 2) mask_y_found = true;
976  }
977 
978  // Found a nearby masked cell (valid n_ind)
979  if (mask_x_found && mask_y_found) {
980  n_ind = std::min(ii,jj) + one;
981  } else if (mask_x_found) {
982  n_ind = ii + one;
983  } else if (mask_y_found) {
984  n_ind = jj + one;
985  // Pesky corner cell
986  } else {
987  if (near_x_lo_wall || near_x_hi_wall) {
988  Real dj_min{width-one};
989  int j_lb = std::max(vbx_lo.y,j-width);
990  int j_ub = std::min(vbx_hi.y,j+width);
991  int li = (near_x_lo_wall) ? vbx_lo.x : vbx_hi.x;
992  for (int lj(j_lb); lj<=j_ub; ++lj) {
993  if (mask_arr(li,lj,k) == 2) {
994  mask_y_found = true;
995  dj_min = std::min(dj_min,(Real) std::abs(lj-j));
996  }
997  }
998  if (mask_y_found) {
999  Real mag = std::sqrt( Real(dj_min*dj_min + ii*ii) );
1000  n_ind = std::min(mag,width-one) + one;
1001  } else {
1002  Abort("Mask not found near x wall!");
1003  }
1004  } else if (near_y_lo_wall || near_y_hi_wall) {
1005  Real di_min{width-one};
1006  int i_lb = std::max(vbx_lo.x,i-width);
1007  int i_ub = std::min(vbx_hi.x,i+width);
1008  int lj = (near_y_lo_wall) ? vbx_lo.y : vbx_hi.y;
1009  for (int li(i_lb); li<=i_ub; ++li) {
1010  if (mask_arr(li,lj,k) == 2) {
1011  mask_x_found = true;
1012  di_min = std::min(di_min,(Real) std::abs(li-i));
1013  }
1014  }
1015  if (mask_x_found) {
1016  Real mag = std::sqrt( Real(di_min*di_min + jj*jj) );
1017  n_ind = std::min(mag,width-one) + one;
1018  } else {
1019  Abort("Mask not found near y wall!");
1020  }
1021  } else {
1022  Abort("Relaxation cell must be near a wall!");
1023  }
1024  }
1025 
1026  Real Factor = (num - n_ind)/denom;
1027  Real d = data_arr(i ,j ,k ,n+icomp) + delta_t*rhs_arr(i , j , k ,n+icomp);
1028  Real d_ip1 = data_arr(i+1,j ,k ,n+icomp) + delta_t*rhs_arr(i+1, j , k ,n+icomp);
1029  Real d_im1 = data_arr(i-1,j ,k ,n+icomp) + delta_t*rhs_arr(i-1, j , k ,n+icomp);
1030  Real d_jp1 = data_arr(i ,j+1,k ,n+icomp) + delta_t*rhs_arr(i , j+1, k ,n+icomp);
1031  Real d_jm1 = data_arr(i ,j-1,k ,n+icomp) + delta_t*rhs_arr(i , j-1, k ,n+icomp);
1032  Real delta = fine_arr(i ,j ,k,n) - d;
1033  Real delta_xp = fine_arr(i+1,j ,k,n) - d_ip1;
1034  Real delta_xm = fine_arr(i-1,j ,k,n) - d_im1;
1035  Real delta_yp = fine_arr(i ,j+1,k,n) - d_jp1;
1036  Real delta_ym = fine_arr(i ,j-1,k,n) - d_jm1;
1037  Real Laplacian = delta_xp + delta_xm + delta_yp + delta_ym - Real(4.0)*delta;
1038  rhs_arr(i,j,k,n) += (F1*delta - F2*Laplacian) * Factor;
1039  }
1040  });
1041  } // mfi
1042  } // ivar_idx
1043 }
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:28
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(fab_box, [=] AMREX_GPU_DEVICE(int i, int j, int k) { qrcuten_arr(i, j, k)=Real(0);qscuten_arr(i, j, k)=Real(0);qicuten_arr(i, j, k)=Real(0);})
VariationOperator::Laplacian
@ Laplacian
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:108
ERFFillPatcher::FillRelax
void FillRelax(amrex::MultiFab &mf, double time, BC &cbc, amrex::Vector< amrex::BCRec > const &bcs)
Definition: ERF_FillPatcher.H:195
ERFFillPatcher::GetSetMaskVal
int GetSetMaskVal()
Definition: ERF_FillPatcher.H:98
ERFFillPatcher::GetRelaxMaskVal
int GetRelaxMaskVal()
Definition: ERF_FillPatcher.H:103
IntVars::NumTypes
@ NumTypes
Definition: ERF_IndexDefines.H:236
IntVars::ymom
@ ymom
Definition: ERF_IndexDefines.H:234
IntVars::cons
@ cons
Definition: ERF_IndexDefines.H:232
IntVars::zmom
@ zmom
Definition: ERF_IndexDefines.H:235
IntVars::xmom
@ xmom
Definition: ERF_IndexDefines.H:233
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)
Here is the call graph for this function:

◆ realbdy_compute_interior_ghost_rhs()

void realbdy_compute_interior_ghost_rhs ( const double &  time,
const double &  delta_t_d,
const double &  start_bdy_time,
const double &  final_bdy_time,
const double &  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,
const bool  use_wrf_bdy_density,
const Real &  bdy_rho_nudge_factor 
)

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
180 {
181  BL_PROFILE_REGION("realbdy_compute_interior_ghost_RHS()");
182 
183  Real delta_t = static_cast<Real>(delta_t_d);
184 
185  //
186  // Note that time (= start_time+old_stage_time) is measured as total time
187  // start_bdy_time and final_bdy_time are also measured as total time
188  //
189 
190  // Get bndry data if we have it
191  Vector<int> bnd_map = {BCVars::xvel_bc, BCVars::yvel_bc, BCVars::RhoTheta_bc_comp};
192  Array4<Real> bdatxlo, bdatxhi, bdatylo, bdatyhi;
193  Array4<Real> btenxlo, btenxhi, btenylo, btenyhi;
194  if (m_r2d) {
195  // Index is [plane orientation] and [level]
196  Vector<std::unique_ptr<PlaneVector>>& bndry_data = m_r2d->interp_in_time(time);
197  bdatxlo = (*bndry_data[0])[0].array();
198  bdatylo = (*bndry_data[1])[0].array();
199  bdatxhi = (*bndry_data[3])[0].array();
200  bdatyhi = (*bndry_data[4])[0].array();
201 
202  Vector<std::unique_ptr<PlaneVector>>& bndry_tend = m_r2d->get_tendency(time);
203  btenxlo = (*bndry_tend[0])[0].array();
204  btenylo = (*bndry_tend[1])[0].array();
205  btenxhi = (*bndry_tend[3])[0].array();
206  btenyhi = (*bndry_tend[4])[0].array();
207  }
208 
209  // Relaxation constants
210  Real F1 = one/(nudge_factor*delta_t);
211  Real F1_rho = one/((bdy_rho_nudge_factor > zero ? bdy_rho_nudge_factor : nudge_factor) * delta_t);
212 
213  // Time interpolation
214  double dT_d = bdy_time_interval;
215 
216  int n_time = static_cast<int>( (time-start_bdy_time) / dT_d);
217  int n_time_p1 = n_time + 1;
218  Real alpha = static_cast<Real>(((time-start_bdy_time) - n_time * dT_d) / dT_d);
219 
220  // Do not over run the last bdy file
221  if (time >= final_bdy_time) {
222  n_time = static_cast<int>( (final_bdy_time - start_bdy_time)/ dT_d);
223  n_time_p1 = n_time;
224  alpha = zero;
225  }
226 
227  AMREX_ALWAYS_ASSERT( alpha >= zero && alpha <= one);
228  Real oma = one - alpha;
229 
230  // Temporary FABs for storage (owned/filled on all ranks)
231  FArrayBox U_xlo, U_xhi, U_ylo, U_yhi;
232  FArrayBox V_xlo, V_xhi, V_ylo, V_yhi;
233  FArrayBox T_xlo, T_xhi, T_ylo, T_yhi;
234  FArrayBox R_xlo, R_xhi, R_ylo, R_yhi;
235  Array4<Real> r_xlo_arr, r_xhi_arr, r_ylo_arr, r_yhi_arr;
236 
237  // Variable index map (WRFBdyVars -> Vars)
238  Vector<int> var_map = {Vars::xvel, Vars::yvel, Vars::cons };
239  Vector<int> ivar_map = {IntVars::xmom, IntVars::ymom, IntVars::cons};
240 
241  // Variable icomp map
242  Vector<int> comp_map = {0, 0, RhoTheta_comp};
243 
244  // Indices
245  int ivarU = RealBdyVars::U;
246  int ivarV = RealBdyVars::V;
247  int ivarT = RealBdyVars::T;
248  // This helper applies the legacy velocity/theta boundary relaxation only.
249  // Density (R) and moisture (QV/QC/QI) are handled by fill_from_realbdy;
250  // do not derive this bound from RealBdyVars::NumTypes because QC/QI are
251  // appended to that serialized layout.
252  int BdyEnd = RealBdyVars::T + 1;
253 
254  // NOTE: The sizing of the temporary BDY FABS is
255  // GLOBAL and occurs over the entire BDY region.
256 
257  // Size the FABs
258  //==========================================================
259  for (int ivar(ivarU); ivar < BdyEnd; ivar++) {
260  int ivar_idx = var_map[ivar];
261  Box domain = geom.Domain();
262  auto ixtype = S_cur_data[ivar_idx].boxArray().ixType();
263  domain.convert(ixtype);
264 
265  // NOTE: Ghost cells needed for idx type mismatch between mask and data (do_upwind)
266  IntVect ng_vect(0);
267  //IntVect ng_vect(1,1,0);
268  Box gdom(domain); gdom.grow(ng_vect);
269  Box bx_xlo, bx_xhi, bx_ylo, bx_yhi;
270  realbdy_interior_bxs_xy(gdom, domain, width,
271  bx_xlo, bx_xhi,
272  bx_ylo, bx_yhi,
273  ng_vect, true);
274 
275  // Size the FABs
276  if (ivar == ivarU) {
277  U_xlo.resize(bx_xlo,1,The_Async_Arena()); U_xhi.resize(bx_xhi,1,The_Async_Arena());
278  U_ylo.resize(bx_ylo,1,The_Async_Arena()); U_yhi.resize(bx_yhi,1,The_Async_Arena());
279  } else if (ivar == ivarV) {
280  V_xlo.resize(bx_xlo,1,The_Async_Arena()); V_xhi.resize(bx_xhi,1,The_Async_Arena());
281  V_ylo.resize(bx_ylo,1,The_Async_Arena()); V_yhi.resize(bx_yhi,1,The_Async_Arena());
282  } else if (ivar == ivarT){
283  T_xlo.resize(bx_xlo,1,The_Async_Arena()); T_xhi.resize(bx_xhi,1,The_Async_Arena());
284  T_ylo.resize(bx_ylo,1,The_Async_Arena()); T_yhi.resize(bx_yhi,1,The_Async_Arena());
285  } else {
286  continue;
287  }
288  } // ivar
289 
290  if (use_wrf_bdy_density) {
291  const IntVect ng_vect(1,1,0);
292  Box domain = geom.Domain();
293  Box gdom(domain); gdom.grow(ng_vect);
294  Box bx_xlo, bx_xhi, bx_ylo, bx_yhi;
295  realbdy_interior_bxs_xy(gdom, domain, width,
296  bx_xlo, bx_xhi, bx_ylo, bx_yhi,
297  ng_vect, true);
298 
299  R_xlo.resize(bx_xlo, 1, The_Async_Arena());
300  R_xhi.resize(bx_xhi, 1, The_Async_Arena());
301  R_ylo.resize(bx_ylo, 1, The_Async_Arena());
302  R_yhi.resize(bx_yhi, 1, The_Async_Arena());
303 
304  const auto& r_xlo_n = bdy_data_xlo[n_time ][WRFBdyVars::R].const_array();
305  const auto& r_xlo_np1 = bdy_data_xlo[n_time_p1][WRFBdyVars::R].const_array();
306  const auto& r_xhi_n = bdy_data_xhi[n_time ][WRFBdyVars::R].const_array();
307  const auto& r_xhi_np1 = bdy_data_xhi[n_time_p1][WRFBdyVars::R].const_array();
308  const auto& r_ylo_n = bdy_data_ylo[n_time ][WRFBdyVars::R].const_array();
309  const auto& r_ylo_np1 = bdy_data_ylo[n_time_p1][WRFBdyVars::R].const_array();
310  const auto& r_yhi_n = bdy_data_yhi[n_time ][WRFBdyVars::R].const_array();
311  const auto& r_yhi_np1 = bdy_data_yhi[n_time_p1][WRFBdyVars::R].const_array();
312  const auto& rbx = lbound(domain);
313  const auto& rhi = ubound(domain);
314  r_xlo_arr = R_xlo.array();
315  r_xhi_arr = R_xhi.array();
316  r_ylo_arr = R_ylo.array();
317  r_yhi_arr = R_yhi.array();
318 
319  ParallelFor(bx_xlo, bx_xhi,
320  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
321  int ii = amrex::max(i, rbx.x); ii = amrex::min(ii, rbx.x + width - 1);
322  int jj = amrex::max(j, rbx.y); jj = amrex::min(jj, rhi.y);
323  r_xlo_arr(i,j,k) = oma*r_xlo_n(ii,jj,k) + alpha*r_xlo_np1(ii,jj,k);
324  },
325  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
326  int ii = amrex::max(i, rhi.x - width + 1); ii = amrex::min(ii, rhi.x);
327  int jj = amrex::max(j, rbx.y); jj = amrex::min(jj, rhi.y);
328  r_xhi_arr(i,j,k) = oma*r_xhi_n(ii,jj,k) + alpha*r_xhi_np1(ii,jj,k);
329  });
330  ParallelFor(bx_ylo, bx_yhi,
331  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
332  int ii = amrex::max(i, rbx.x); ii = amrex::min(ii, rhi.x);
333  int jj = amrex::max(j, rbx.y); jj = amrex::min(jj, rbx.y + width - 1);
334  r_ylo_arr(i,j,k) = oma*r_ylo_n(ii,jj,k) + alpha*r_ylo_np1(ii,jj,k);
335  },
336  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
337  int ii = amrex::max(i, rbx.x); ii = amrex::min(ii, rhi.x);
338  int jj = amrex::max(j, rhi.y - width + 1); jj = amrex::min(jj, rhi.y);
339  r_yhi_arr(i,j,k) = oma*r_yhi_n(ii,jj,k) + alpha*r_yhi_np1(ii,jj,k);
340  });
341  }
342 
343 
344  // NOTE: These operations use the BDY FABS and RHO. The
345  // use of RHO to go from PRIM -> CONS requires that
346  // these operations be LOCAL. So we have allocated
347  // enough space to do global operations (1 rank) but
348  // will fill a subset of that data that the rank owns.
349 
350  // Populate FABs from bdy interpolation (primitive vars)
351  //==========================================================
352  for (int ivar(ivarU); ivar < BdyEnd; ivar++) {
353  int ivar_idx = var_map[ivar];
354  Box domain = geom.Domain();
355  auto ixtype = S_cur_data[ivar_idx].boxArray().ixType();
356  domain.convert(ixtype);
357  const auto& dom_lo = lbound(domain);
358  const auto& dom_hi = ubound(domain);
359 
360  // BndryReg idx and limiting
361  int bdy_comp = bnd_map[ivar];
362  const auto& dom_cc_lo = lbound(geom.Domain());
363  const auto& dom_cc_hi = ubound(geom.Domain());
364 
365 #ifdef _OPENMP
366 #pragma omp parallel if (Gpu::notInLaunchRegion())
367 #endif
368  for (MFIter mfi(S_cur_data[ivar_idx],TilingIfNotGPU()); mfi.isValid(); ++mfi) {
369  // NOTE: Ghost cells needed for idx type mismatch between mask and data (do_upwind)
370  IntVect ng_vect(0);
371  //IntVect ng_vect(1,1,0);
372  Box gtbx = grow(mfi.tilebox(ixtype.toIntVect()),ng_vect);
373  Box tbx_xlo, tbx_xhi, tbx_ylo, tbx_yhi;
374  realbdy_interior_bxs_xy(gtbx, domain, width,
375  tbx_xlo, tbx_xhi,
376  tbx_ylo, tbx_yhi,
377  ng_vect, true);
378 
379  Array4<Real> arr_xlo; Array4<Real> arr_xhi;
380  Array4<Real> arr_ylo; Array4<Real> arr_yhi;
381  if (ivar == ivarU) {
382  arr_xlo = U_xlo.array(); arr_xhi = U_xhi.array();
383  arr_ylo = U_ylo.array(); arr_yhi = U_yhi.array();
384  } else if (ivar == ivarV) {
385  arr_xlo = V_xlo.array(); arr_xhi = V_xhi.array();
386  arr_ylo = V_ylo.array(); arr_yhi = V_yhi.array();
387  } else if (ivar == ivarT){
388  arr_xlo = T_xlo.array(); arr_xhi = T_xhi.array();
389  arr_ylo = T_ylo.array(); arr_yhi = T_yhi.array();
390  } else {
391  continue;
392  }
393 
394  // Boundary data at fixed time intervals
395  const auto& bdatxlo_n = bdy_data_xlo[n_time ][ivar].const_array();
396  const auto& bdatxlo_np1 = bdy_data_xlo[n_time_p1][ivar].const_array();
397  const auto& bdatxhi_n = bdy_data_xhi[n_time ][ivar].const_array();
398  const auto& bdatxhi_np1 = bdy_data_xhi[n_time_p1][ivar].const_array();
399  const auto& bdatylo_n = bdy_data_ylo[n_time ][ivar].const_array();
400  const auto& bdatylo_np1 = bdy_data_ylo[n_time_p1][ivar].const_array();
401  const auto& bdatyhi_n = bdy_data_yhi[n_time ][ivar].const_array();
402  const auto& bdatyhi_np1 = bdy_data_yhi[n_time_p1][ivar].const_array();
403 
404  // Current density to convert to conserved vars
405  Array4<Real> r_arr = S_cur_data[IntVars::cons].array(mfi);
406 
407  // Limiting offset
408  int offset = width - 1;
409 
410  // Populate with interpolation (protect from ghost cells)
411  const auto rxlo = r_xlo_arr;
412  const auto rxhi = r_xhi_arr;
413  const auto rylo = r_ylo_arr;
414  const auto ryhi = r_yhi_arr;
415  ParallelFor(tbx_xlo, tbx_xhi,
416  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
417  {
418  int ii = std::max(i , dom_lo.x); ii = std::min(ii, dom_lo.x+offset);
419  int jj = std::max(j , dom_lo.y); jj = std::min(jj, dom_hi.y);
420 
421  Real rho_interp;
422  if (use_wrf_bdy_density && ivar==ivarU) {
423  int im = amrex::max(i-1, dom_lo.x);
424  rho_interp = myhalf * (rxlo(im,j,k) + rxlo(amrex::max(i, dom_lo.x),j,k));
425  } else if (use_wrf_bdy_density && ivar==ivarV) {
426  int jm = amrex::max(j-1, dom_lo.y);
427  rho_interp = myhalf * (rxlo(i,jm,k) + rxlo(i,amrex::max(j, dom_lo.y),k));
428  } else if (use_wrf_bdy_density) {
429  rho_interp = rxlo(i,j,k);
430  } else if (ivar==ivarU) {
431  rho_interp = myhalf * ( r_arr(i-1,j ,k) + r_arr(i,j,k) );
432  } else if (ivar==ivarV) {
433  rho_interp = myhalf * ( r_arr(i ,j-1,k) + r_arr(i,j,k) );
434  } else {
435  rho_interp = r_arr(i,j,k);
436  }
437 
438  if (bdatxlo) {
439  int ii2 = std::min(std::max(i , dom_cc_lo.x), dom_cc_hi.x);
440  int jj2 = std::min(std::max(j , dom_cc_lo.y), dom_cc_hi.y);
441  arr_xlo(i,j,k) = rho_interp * bdatxlo(ii2,jj2,k,bdy_comp);
442  } else {
443  arr_xlo(i,j,k) = rho_interp * ( oma * bdatxlo_n (ii,jj,k,0)
444  + alpha * bdatxlo_np1(ii,jj,k,0) );
445  }
446  },
447  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
448  {
449  int ii = std::max(i , dom_hi.x-offset); ii = std::min(ii, dom_hi.x);
450  int jj = std::max(j , dom_lo.y); jj = std::min(jj, dom_hi.y);
451 
452  Real rho_interp;
453  if (use_wrf_bdy_density && ivar==ivarU) {
454  int im = amrex::max(i-1, dom_lo.x);
455  rho_interp = myhalf * (rxhi(im,j,k) + rxhi(amrex::max(i, dom_lo.x),j,k));
456  } else if (use_wrf_bdy_density && ivar==ivarV) {
457  int jm = amrex::max(j-1, dom_lo.y);
458  rho_interp = myhalf * (rxhi(i,jm,k) + rxhi(i,amrex::max(j, dom_lo.y),k));
459  } else if (use_wrf_bdy_density) {
460  rho_interp = rxhi(i,j,k);
461  } else if (ivar==ivarU) {
462  rho_interp = myhalf * ( r_arr(i-1,j ,k) + r_arr(i,j,k) );
463  } else if (ivar==ivarV) {
464  rho_interp = myhalf * ( r_arr(i ,j-1,k) + r_arr(i,j,k) );
465  } else {
466  rho_interp = r_arr(i,j,k);
467  }
468 
469  if (bdatxhi) {
470  int ii2 = std::min(std::max(i , dom_cc_lo.x), dom_cc_hi.x);
471  int jj2 = std::min(std::max(j , dom_cc_lo.y), dom_cc_hi.y);
472  arr_xhi(i,j,k) = rho_interp * bdatxhi(ii2,jj2,k,bdy_comp);
473  } else {
474  arr_xhi(i,j,k) = rho_interp * ( oma * bdatxhi_n (ii,jj,k,0)
475  + alpha * bdatxhi_np1(ii,jj,k,0) );
476  }
477  });
478 
479  ParallelFor(tbx_ylo, tbx_yhi,
480  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
481  {
482  int ii = std::max(i , dom_lo.x); ii = std::min(ii, dom_hi.x);
483  int jj = std::max(j , dom_lo.y); jj = std::min(jj, dom_lo.y+offset);
484 
485  Real rho_interp;
486  if (use_wrf_bdy_density && ivar==ivarU) {
487  int im = amrex::max(i-1, dom_lo.x);
488  rho_interp = myhalf * (rylo(im,j,k) + rylo(amrex::max(i, dom_lo.x),j,k));
489  } else if (use_wrf_bdy_density && ivar==ivarV) {
490  int jm = amrex::max(j-1, dom_lo.y);
491  rho_interp = myhalf * (rylo(i,jm,k) + rylo(i,amrex::max(j, dom_lo.y),k));
492  } else if (use_wrf_bdy_density) {
493  rho_interp = rylo(i,j,k);
494  } else if (ivar==ivarU) {
495  rho_interp = myhalf * ( r_arr(i-1,j ,k) + r_arr(i,j,k) );
496  } else if (ivar==ivarV) {
497  rho_interp = myhalf * ( r_arr(i ,j-1,k) + r_arr(i,j,k) );
498  } else {
499  rho_interp = r_arr(i,j,k);
500  }
501 
502  if (bdatylo) {
503  int ii2 = std::min(std::max(i , dom_cc_lo.x), dom_cc_hi.x);
504  int jj2 = std::min(std::max(j , dom_cc_lo.y), dom_cc_hi.y);
505  arr_ylo(i,j,k) = rho_interp * bdatylo(ii2,jj2,k,bdy_comp);
506  } else {
507  arr_ylo(i,j,k) = rho_interp * ( oma * bdatylo_n (ii,jj,k,0)
508  + alpha * bdatylo_np1(ii,jj,k,0) );
509  }
510  },
511  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
512  {
513  int ii = std::max(i , dom_lo.x); ii = std::min(ii, dom_hi.x);
514  int jj = std::max(j , dom_hi.y-offset); jj = std::min(jj, dom_hi.y);
515 
516  Real rho_interp;
517  if (use_wrf_bdy_density && ivar==ivarU) {
518  int im = amrex::max(i-1, dom_lo.x);
519  rho_interp = myhalf * (ryhi(im,j,k) + ryhi(amrex::max(i, dom_lo.x),j,k));
520  } else if (use_wrf_bdy_density && ivar==ivarV) {
521  int jm = amrex::max(j-1, dom_lo.y);
522  rho_interp = myhalf * (ryhi(i,jm,k) + ryhi(i,amrex::max(j, dom_lo.y),k));
523  } else if (use_wrf_bdy_density) {
524  rho_interp = ryhi(i,j,k);
525  } else if (ivar==ivarU) {
526  rho_interp = myhalf * ( r_arr(i-1,j ,k) + r_arr(i,j,k) );
527  } else if (ivar==ivarV) {
528  rho_interp = myhalf * ( r_arr(i ,j-1,k) + r_arr(i,j,k) );
529  } else {
530  rho_interp = r_arr(i,j,k);
531  }
532 
533  if (bdatyhi) {
534  int ii2 = std::min(std::max(i , dom_cc_lo.x), dom_cc_hi.x);
535  int jj2 = std::min(std::max(j , dom_cc_lo.y), dom_cc_hi.y);
536  arr_yhi(i,j,k) = rho_interp * bdatyhi(ii2,jj2,k,bdy_comp);
537  } else {
538  arr_yhi(i,j,k) = rho_interp * ( oma * bdatyhi_n (ii,jj,k,0)
539  + alpha * bdatyhi_np1(ii,jj,k,0) );
540  }
541  });
542  } // mfi
543  } // ivar
544 
545 
546  // Compute RHS in relaxation region
547  //==========================================================
548  auto dx = geom.CellSizeArray();
549  auto ProbLo = geom.ProbLoArray();
550  auto ProbHi = geom.ProbHiArray();
551 
552  for (int ivar(ivarU); ivar < BdyEnd; ivar++) {
553  int ivar_idx = ivar_map[ivar];
554  int icomp = comp_map[ivar];
555 
556  Box domain = geom.Domain();
557  domain.convert(S_cur_data[ivar_idx].boxArray().ixType());
558  IntVect ng_vect(0);
559 
560 #ifdef _OPENMP
561 #pragma omp parallel if (Gpu::notInLaunchRegion())
562 #endif
563  for (MFIter mfi(S_cur_data[ivar_idx],TilingIfNotGPU()); mfi.isValid(); ++mfi)
564  {
565  Box tbx = mfi.tilebox();
566  Box tbx_xlo, tbx_xhi, tbx_ylo, tbx_yhi;
567  realbdy_interior_bxs_xy(tbx, domain, width,
568  tbx_xlo, tbx_xhi,
569  tbx_ylo, tbx_yhi,
570  ng_vect);
571 
572  Array4<Real> rhs_arr; Array4<Real> data_arr;
573  Array4<Real> arr_xlo; Array4<Real> arr_xhi;
574  Array4<Real> arr_ylo; Array4<Real> arr_yhi;
575  if (ivar == ivarU) {
576  arr_xlo = U_xlo.array(); arr_xhi = U_xhi.array();
577  arr_ylo = U_ylo.array(); arr_yhi = U_yhi.array();
578  rhs_arr = S_rhs[IntVars::xmom].array(mfi);
579  data_arr = S_cur_data[IntVars::xmom].array(mfi);
580  } else if (ivar == ivarV) {
581  arr_xlo = V_xlo.array(); arr_xhi = V_xhi.array();
582  arr_ylo = V_ylo.array(); arr_yhi = V_yhi.array();
583  rhs_arr = S_rhs[IntVars::ymom].array(mfi);
584  data_arr = S_cur_data[IntVars::ymom].array(mfi);
585  } else if (ivar == ivarT){
586  arr_xlo = T_xlo.array(); arr_xhi = T_xhi.array();
587  arr_ylo = T_ylo.array(); arr_yhi = T_yhi.array();
588  rhs_arr = S_rhs[IntVars::cons].array(mfi);
589  data_arr = S_cur_data[IntVars::cons].array(mfi);
590  } else {
591  continue;
592  }
593 
594  realbdy_compute_relaxation(icomp, 1,
595  width, dx, ProbLo, ProbHi, F1,
596  tbx_xlo , tbx_xhi , tbx_ylo , tbx_yhi ,
597  arr_xlo , arr_xhi , arr_ylo , arr_yhi ,
598  data_arr, rhs_arr , c_p, rdOcp);
599  } // mfi
600  } // ivar
601 
602  if (use_wrf_bdy_density) {
603  for (MFIter mfi(S_cur_data[IntVars::cons], TilingIfNotGPU()); mfi.isValid(); ++mfi) {
604  Box domain = geom.Domain();
605  Box tbx_xlo, tbx_xhi, tbx_ylo, tbx_yhi;
606  realbdy_interior_bxs_xy(mfi.tilebox(), domain, width,
607  tbx_xlo, tbx_xhi, tbx_ylo, tbx_yhi);
608  realbdy_compute_relaxation(Rho_comp, 1, width, dx, ProbLo, ProbHi, F1_rho,
609  tbx_xlo, tbx_xhi, tbx_ylo, tbx_yhi,
610  r_xlo_arr, r_xhi_arr, r_ylo_arr, r_yhi_arr,
611  S_cur_data[IntVars::cons].const_array(mfi),
612  S_rhs[IntVars::cons].array(mfi), c_p, rdOcp);
613  }
614  }
615 
616  // Set normal velocity RHS at the boundary
617  //==========================================================
618  Box domain = geom.Domain();
619  Box domainx = convert(domain, IntVect(1,0,0));
620  Box domainy = convert(domain, IntVect(0,1,0));
621 
622  int ilo = domainx.smallEnd(0);
623  int ihi = domainx.bigEnd(0);
624  int jlo = domainy.smallEnd(1);
625  int jhi = domainy.bigEnd(1);
626 
627 #ifdef _OPENMP
628 #pragma omp parallel if (Gpu::notInLaunchRegion())
629 #endif
630  for (MFIter mfi(S_cur_data[IntVars::cons],TilingIfNotGPU()); mfi.isValid(); ++mfi)
631  {
632  Box tbx = mfi.nodaltilebox(0);
633  Box tbx_lo, tbx_hi;
634  if (tbx.smallEnd(0) == ilo) {
635  tbx_lo = makeSlab(tbx,0,ilo);
636  }
637  if (tbx.bigEnd(0) == ihi) {
638  tbx_hi = makeSlab(tbx,0,ihi);
639  }
640 
641  Box tby = mfi.nodaltilebox(1);
642  Box tby_lo, tby_hi;
643  if (tby.smallEnd(1) == jlo) {
644  tby_lo = makeSlab(tby,1,jlo);
645  }
646  if (tby.bigEnd(1) == jhi) {
647  tby_hi = makeSlab(tby,1,jhi);
648  }
649 
650  Array4<Real> rhs_xmom = S_rhs[IntVars::xmom].array(mfi);
651  Array4<Real> rhs_ymom = S_rhs[IntVars::ymom].array(mfi);
652 
653  Array4<const Real> rhs_cons = S_rhs[IntVars::cons].const_array(mfi);
654  Array4<const Real> cons_arr = S_cur_data[IntVars::cons].const_array(mfi);
655 
656  const auto& bdatxlo_n = bdy_data_xlo[n_time ][ivarU].const_array();
657  const auto& bdatxlo_np1 = bdy_data_xlo[n_time_p1][ivarU].const_array();
658  const auto& bdatxhi_n = bdy_data_xhi[n_time ][ivarU].const_array();
659  const auto& bdatxhi_np1 = bdy_data_xhi[n_time_p1][ivarU].const_array();
660 
661  const auto& bdatylo_n = bdy_data_ylo[n_time ][ivarV].const_array();
662  const auto& bdatylo_np1 = bdy_data_ylo[n_time_p1][ivarV].const_array();
663  const auto& bdatyhi_n = bdy_data_yhi[n_time ][ivarV].const_array();
664  const auto& bdatyhi_np1 = bdy_data_yhi[n_time_p1][ivarV].const_array();
665 
666  Real dT = static_cast<Real>(bdy_time_interval);
667 
668  ParallelFor(tbx_lo, tbx_hi,
669  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
670  {
671  Real rho_tend = rhs_cons(i,j,k);
672  Real rho_val = Real(0.5) * (cons_arr(i,j,k) + cons_arr(i-1,j,k));
673  Real u_tend, u_val;
674  if (btenxlo) {
675  u_tend = btenxlo(i,j,k,BCVars::xvel_bc);
676  u_val = bdatxlo(i,j,k,BCVars::xvel_bc);
677  } else {
678  u_tend = (bdatxlo_np1(i,j,k) - bdatxlo_n(i,j,k)) / dT;
679  u_val = oma * bdatxlo_n(i,j,k) + alpha * bdatxlo_np1(i,j,k);
680  }
681  rhs_xmom(i,j,k) = rho_val * u_tend + u_val * rho_tend;
682  },
683  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
684  {
685  // NOTE: (i,j,k) is the hi face of the domain, so the cell (i,j,k) is
686  // an exterior ghost cell for which the cons RHS is never filled;
687  // use the adjacent interior cell, mirroring the lo side.
688  Real rho_tend = rhs_cons(i-1,j,k);
689  Real rho_val = Real(0.5) * (cons_arr(i,j,k) + cons_arr(i-1,j,k));
690  Real u_tend, u_val;
691  if (btenxhi) {
692  u_tend = btenxhi(i,j,k,BCVars::xvel_bc);
693  u_val = bdatxhi(i,j,k,BCVars::xvel_bc);
694  } else {
695  u_tend = (bdatxhi_np1(i,j,k) - bdatxhi_n(i,j,k)) / dT;
696  u_val = oma * bdatxhi_n(i,j,k) + alpha * bdatxhi_np1(i,j,k);
697  }
698  rhs_xmom(i,j,k) = rho_val * u_tend + u_val * rho_tend;
699  });
700 
701  ParallelFor(tby_lo, tby_hi,
702  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
703  {
704  Real rho_tend = rhs_cons(i,j,k);
705  Real rho_val = Real(0.5) * (cons_arr(i,j,k) + cons_arr(i,j-1,k));
706  Real v_tend, v_val;
707  if (btenylo) {
708  v_tend = btenylo(i,j,k,BCVars::yvel_bc);
709  v_val = bdatylo(i,j,k,BCVars::yvel_bc);
710  } else {
711  v_tend = (bdatylo_np1(i,j,k) - bdatylo_n(i,j,k)) / dT;
712  v_val = oma * bdatylo_n(i,j,k) + alpha * bdatylo_np1(i,j,k);
713  }
714  rhs_ymom(i,j,k) = rho_val * v_tend + v_val * rho_tend;
715  },
716  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
717  {
718  // NOTE: (i,j,k) is the hi face of the domain, so the cell (i,j,k) is
719  // an exterior ghost cell for which the cons RHS is never filled;
720  // use the adjacent interior cell, mirroring the lo side.
721  Real rho_tend = rhs_cons(i,j-1,k);
722  Real rho_val = Real(0.5) * (cons_arr(i,j,k) + cons_arr(i,j-1,k));
723  Real v_tend, v_val;
724  if (btenyhi) {
725  v_tend = btenyhi(i,j,k,BCVars::yvel_bc);
726  v_val = bdatyhi(i,j,k,BCVars::yvel_bc);
727  } else {
728  v_tend = (bdatyhi_np1(i,j,k) - bdatyhi_n(i,j,k)) / dT;
729  v_val = oma * bdatyhi_n(i,j,k) + alpha * bdatyhi_np1(i,j,k);
730  }
731  rhs_ymom(i,j,k) = rho_val * v_tend + v_val * rho_tend;
732  });
733  } // mfi
734 }
Rho_comp
#define Rho_comp
Definition: ERF_IndexDefines.H:39
RhoTheta_comp
#define RhoTheta_comp
Definition: ERF_IndexDefines.H:40
rdOcp
const Real rdOcp
Definition: ERF_InitCustomPert_ABL.H:72
dx
const Real dx
Definition: ERF_InitCustomPert_ABL.H:44
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:31
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:329
BCVars::RhoTheta_bc_comp
@ RhoTheta_bc_comp
Definition: ERF_IndexDefines.H:91
BCVars::yvel_bc
@ yvel_bc
Definition: ERF_IndexDefines.H:106
BCVars::xvel_bc
@ xvel_bc
Definition: ERF_IndexDefines.H:105
RealBdyVars::U
@ U
Definition: ERF_IndexDefines.H:126
RealBdyVars::T
@ T
Definition: ERF_IndexDefines.H:128
RealBdyVars::V
@ V
Definition: ERF_IndexDefines.H:127
Vars::xvel
@ xvel
Definition: ERF_IndexDefines.H:215
Vars::cons
@ cons
Definition: ERF_IndexDefines.H:214
Vars::yvel
@ yvel
Definition: ERF_IndexDefines.H:216
WRFBdyVars::R
@ R
Definition: ERF_IndexDefines.H:155
WSM6Ind::rhi
@ rhi
Definition: ERF_WSM6.H:250
mp_wdm6::alpha
real(kind=kind_phys), parameter, private alpha
Definition: ERF_module_mp_wdm6.F90:62
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().