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 compute_interior_ghost_bxs_xy (const Box &bx, const Box &domain, const int &width, const int &set_width, Box &bx_xlo, Box &bx_xhi, Box &bx_ylo, Box &bx_yhi, const IntVect &ng_vect, const bool get_int_ng)
 
void realbdy_compute_interior_ghost_rhs (const Real &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

◆ compute_interior_ghost_bxs_xy()

void compute_interior_ghost_bxs_xy ( const Box &  bx,
const Box &  domain,
const int &  width,
const int &  set_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
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  // The four boxes surrounding 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  //==================================================================
59  // NOTE: 4 boxes now encompass interior cells with thickness 'width'
60  //==================================================================
61 
62  gdom_xlo.growLo(0,-set_width); gdom_xhi.growHi(0,-set_width);
63  gdom_xlo.grow (1,-set_width); gdom_xhi.grow (1,-set_width);
64  gdom_ylo.growLo(1,-set_width); gdom_yhi.growHi(1,-set_width);
65 
66  //==================================================================
67  // NOTE: 4 boxes now exclude the regions being set by bndry
68  //==================================================================
69 
70  // Grow boxes to get external ghost cells only
71  gdom_xlo.growLo(0,ng_vect[0]); gdom_xhi.growHi(0,ng_vect[0]);
72  gdom_xlo.grow (1,ng_vect[1]); gdom_xhi.grow (1,ng_vect[1]);
73  gdom_ylo.growLo(1,ng_vect[1]); gdom_yhi.growHi(1,ng_vect[1]);
74 
75  // Grow boxes to get internal ghost cells
76  if (get_int_ng) {
77  gdom_xlo.growHi(0,ng_vect[0]); gdom_xhi.growLo(0,ng_vect[0]);
78  gdom_ylo.grow (0,ng_vect[0]); gdom_yhi.grow (0,ng_vect[0]);
79  gdom_ylo.growHi(1,ng_vect[1]); gdom_yhi.growLo(1,ng_vect[1]);
80  }
81 
82  // Populate everything
83  bx_xlo = (bx & gdom_xlo);
84  bx_xhi = (bx & gdom_xhi);
85  bx_ylo = (bx & gdom_ylo);
86  bx_yhi = (bx & gdom_yhi);
87 }

Referenced by realbdy_compute_interior_ghost_rhs().

Here is the caller graph for this function:

◆ 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
532 {
533  BL_PROFILE_REGION("fine_compute_interior_ghost_RHS()");
534 
535  // Relaxation constants
536  Real F1 = 1./(10.*delta_t);
537  Real F2 = 1./(50.*delta_t);
538 
539  // Vector of MFs to hold data (dm differs w/ fine patch)
540  Vector<MultiFab> fmf_p_v;
541 
542  // Loop over the variables
543  for (int ivar_idx = 0; ivar_idx < IntVars::NumTypes; ++ivar_idx)
544  {
545  // Fine mfs
546  MultiFab& fmf = S_data_f[ivar_idx];
547  MultiFab& rhs = S_rhs_f [ivar_idx];
548 
549  // NOTE: These temporary MFs and copy operations are horrible
550  // for memory usage and efficiency. However, we need to
551  // have access to ghost cells in the cons array to convert
552  // from primitive u/v/w to momentum. Furthermore, the BA
553  // for the fine patches in ERFFillPatcher don't match the
554  // BA for the data/RHS. For this reason, the data is copied
555  // to a vector of MFs (with ghost cells) so the BAs match
556  // the BA of data/RHS and we have access to rho to convert
557  // prim to conserved.
558 
559  // Temp MF on box (distribution map differs w/ fine patch)
560  int num_var = fmf.nComp();
561  fmf_p_v.emplace_back(fmf.boxArray(), fmf.DistributionMap(), num_var, fmf.nGrowVect());
562  MultiFab& fmf_p = fmf_p_v[ivar_idx];
563  MultiFab::Copy(fmf_p,fmf, 0, 0, num_var, fmf.nGrowVect());
564 
565  // Integer mask MF
566  int set_mask_val;
567  int relax_mask_val;
568  iMultiFab* mask;
569 
570  // Fill fine patch on interior halo region
571  //==========================================================
572  if (ivar_idx == IntVars::cons)
573  {
574  FPr_c->FillRelax(fmf_p, time, void_bc, domain_bcs_type);
575  mask = FPr_c->GetMask();
576  set_mask_val = FPr_c->GetSetMaskVal();
577  relax_mask_val = FPr_c->GetRelaxMaskVal();
578  }
579  else if (ivar_idx == IntVars::xmom)
580  {
581  FPr_u->FillRelax(fmf_p, time, void_bc, domain_bcs_type);
582  mask = FPr_u->GetMask();
583  set_mask_val = FPr_u->GetSetMaskVal();
584  relax_mask_val = FPr_u->GetRelaxMaskVal();
585 
586 #ifdef _OPENMP
587 #pragma omp parallel if (Gpu::notInLaunchRegion())
588 #endif
589  for ( MFIter mfi(fmf_p,TilingIfNotGPU()); mfi.isValid(); ++mfi)
590  {
591  Box vbx = mfi.validbox();
592 
593  const Array4<Real>& prim_arr = fmf_p.array(mfi);
594  const Array4<const Real>& rho_arr = fmf_p_v[0].const_array(mfi);
595  const Array4<const int>& mask_arr = mask->const_array(mfi);
596 
597  ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
598  {
599  if (mask_arr(i,j,k) == relax_mask_val) {
600  Real rho_interp = 0.5 * ( rho_arr(i-1,j,k) + rho_arr(i,j,k) );
601  prim_arr(i,j,k) *= rho_interp;
602  }
603  });
604  } // mfi
605  }
606  else if (ivar_idx == IntVars::ymom)
607  {
608  FPr_v->FillRelax(fmf_p, time, void_bc, domain_bcs_type);
609  mask = FPr_v->GetMask();
610  set_mask_val = FPr_v->GetSetMaskVal();
611  relax_mask_val = FPr_v->GetRelaxMaskVal();
612 
613 #ifdef _OPENMP
614 #pragma omp parallel if (Gpu::notInLaunchRegion())
615 #endif
616  for ( MFIter mfi(fmf_p,TilingIfNotGPU()); mfi.isValid(); ++mfi)
617  {
618  Box vbx = mfi.validbox();
619 
620  const Array4<Real>& prim_arr = fmf_p.array(mfi);
621  const Array4<const Real>& rho_arr = fmf_p_v[0].const_array(mfi);
622  const Array4<const int>& mask_arr = mask->const_array(mfi);
623 
624  ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
625  {
626  if (mask_arr(i,j,k) == relax_mask_val) {
627  Real rho_interp = 0.5 * ( rho_arr(i,j-1,k) + rho_arr(i,j,k) );
628  prim_arr(i,j,k) *= rho_interp;
629  }
630  });
631  } // mfi
632  }
633  else if (ivar_idx == IntVars::zmom)
634  {
635  FPr_w->FillRelax(fmf_p, time, void_bc, domain_bcs_type);
636  mask = FPr_w->GetMask();
637  set_mask_val = FPr_w->GetSetMaskVal();
638  relax_mask_val = FPr_w->GetRelaxMaskVal();
639 
640 #ifdef _OPENMP
641 #pragma omp parallel if (Gpu::notInLaunchRegion())
642 #endif
643  for ( MFIter mfi(fmf_p,TilingIfNotGPU()); mfi.isValid(); ++mfi)
644  {
645  Box vbx = mfi.validbox();
646 
647  const Array4<Real>& prim_arr = fmf_p.array(mfi);
648  const Array4<const Real>& rho_arr = fmf_p_v[0].const_array(mfi);
649  const Array4<const int>& mask_arr = mask->const_array(mfi);
650 
651  ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
652  {
653  if (mask_arr(i,j,k) == relax_mask_val) {
654  Real rho_interp = 0.5 * ( rho_arr(i,j,k-1) + rho_arr(i,j,k) );
655  prim_arr(i,j,k) *= rho_interp;
656  }
657  });
658  } // mfi
659  } else {
660  Abort("Dont recognize this variable type in fine_compute_interior_ghost_RHS");
661  }
662 
663 
664  // Zero RHS in set region
665  //==========================================================
666 #ifdef _OPENMP
667 #pragma omp parallel if (Gpu::notInLaunchRegion())
668 #endif
669  for ( MFIter mfi(rhs,TilingIfNotGPU()); mfi.isValid(); ++mfi)
670  {
671  Box vbx = mfi.validbox();
672  const Array4<Real>& rhs_arr = rhs.array(mfi);
673  const Array4<const int>& mask_arr = mask->const_array(mfi);
674 
675  ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
676  {
677  if (mask_arr(i,j,k) == set_mask_val) {
678  rhs_arr(i,j,k) = 0.0;
679  }
680  });
681  } // mfi
682 
683  // For Laplacian stencil
684  rhs.FillBoundary(geom.periodicity());
685 
686 
687  // Compute RHS in relaxation region
688  //==========================================================
689 #ifdef _OPENMP
690 #pragma omp parallel if (Gpu::notInLaunchRegion())
691 #endif
692  for ( MFIter mfi(fmf_p,TilingIfNotGPU()); mfi.isValid(); ++mfi)
693  {
694  Box vbx = mfi.validbox();
695  const Array4<Real>& rhs_arr = rhs.array(mfi);
696  const Array4<const Real>& fine_arr = fmf_p.const_array(mfi);
697  const Array4<const Real>& data_arr = fmf.const_array(mfi);
698  const Array4<const int>& mask_arr = mask->const_array(mfi);
699 
700  const auto& vbx_lo = lbound(vbx);
701  const auto& vbx_hi = ubound(vbx);
702 
703  int icomp = 0;
704 
705  int Spec_z = set_width;
706  int Relax_z = width - Spec_z;
707  Real num = Real(Spec_z + Relax_z);
708  Real denom = Real(Relax_z - 1);
709  ParallelFor(vbx, num_var, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
710  {
711  if (mask_arr(i,j,k) == relax_mask_val) {
712 
713  // Indices
714  Real n_ind(-1); // Set to -1 to quiet compiler warning
715  int ii(width-1); int jj(width-1);
716  bool near_x_lo_wall(false); bool near_x_hi_wall(false);
717  bool near_y_lo_wall(false); bool near_y_hi_wall(false);
718  bool mask_x_found(false); bool mask_y_found(false);
719 
720  // Near x-wall
721  if ((i-vbx_lo.x) < width) {
722  near_x_lo_wall = true;
723  ii = i-vbx_lo.x;
724  if (mask_arr(vbx_lo.x,j,k) == 2) mask_x_found = true;
725  } else if ((vbx_hi.x-i) < width) {
726  near_x_hi_wall = true;
727  ii = vbx_hi.x-i;
728  if (mask_arr(vbx_hi.x,j,k) == 2) mask_x_found = true;
729  }
730 
731  // Near y-wall
732  if ((j-vbx_lo.y) < width) {
733  near_y_lo_wall = true;
734  jj = j-vbx_lo.y;
735  if (mask_arr(i,vbx_lo.y,k) == 2) mask_y_found = true;
736  } else if ((vbx_hi.y-j) < width) {
737  near_y_hi_wall = true;
738  jj = vbx_hi.y-j;
739  if (mask_arr(i,vbx_hi.y,k) == 2) mask_y_found = true;
740  }
741 
742  // Found a nearby masked cell (valid n_ind)
743  if (mask_x_found && mask_y_found) {
744  n_ind = std::min(ii,jj) + 1.0;
745  } else if (mask_x_found) {
746  n_ind = ii + 1.0;
747  } else if (mask_y_found) {
748  n_ind = jj + 1.0;
749  // Pesky corner cell
750  } else {
751  if (near_x_lo_wall || near_x_hi_wall) {
752  Real dj_min{width-1.0};
753  int j_lb = std::max(vbx_lo.y,j-width);
754  int j_ub = std::min(vbx_hi.y,j+width);
755  int li = (near_x_lo_wall) ? vbx_lo.x : vbx_hi.x;
756  for (int lj(j_lb); lj<=j_ub; ++lj) {
757  if (mask_arr(li,lj,k) == 2) {
758  mask_y_found = true;
759  dj_min = std::min(dj_min,(Real) std::abs(lj-j));
760  }
761  }
762  if (mask_y_found) {
763  Real mag = sqrt( Real(dj_min*dj_min + ii*ii) );
764  n_ind = std::min(mag,width-1.0) + 1.0;
765  } else {
766  Abort("Mask not found near x wall!");
767  }
768  } else if (near_y_lo_wall || near_y_hi_wall) {
769  Real di_min{width-1.0};
770  int i_lb = std::max(vbx_lo.x,i-width);
771  int i_ub = std::min(vbx_hi.x,i+width);
772  int lj = (near_y_lo_wall) ? vbx_lo.y : vbx_hi.y;
773  for (int li(i_lb); li<=i_ub; ++li) {
774  if (mask_arr(li,lj,k) == 2) {
775  mask_x_found = true;
776  di_min = std::min(di_min,(Real) std::abs(li-i));
777  }
778  }
779  if (mask_x_found) {
780  Real mag = sqrt( Real(di_min*di_min + jj*jj) );
781  n_ind = std::min(mag,width-1.0) + 1.0;
782  } else {
783  Abort("Mask not found near y wall!");
784  }
785  } else {
786  Abort("Relaxation cell must be near a wall!");
787  }
788  }
789 
790  Real Factor = (num - n_ind)/denom;
791  Real d = data_arr(i ,j ,k ,n+icomp) + delta_t*rhs_arr(i , j , k ,n+icomp);
792  Real d_ip1 = data_arr(i+1,j ,k ,n+icomp) + delta_t*rhs_arr(i+1, j , k ,n+icomp);
793  Real d_im1 = data_arr(i-1,j ,k ,n+icomp) + delta_t*rhs_arr(i-1, j , k ,n+icomp);
794  Real d_jp1 = data_arr(i ,j+1,k ,n+icomp) + delta_t*rhs_arr(i , j+1, k ,n+icomp);
795  Real d_jm1 = data_arr(i ,j-1,k ,n+icomp) + delta_t*rhs_arr(i , j-1, k ,n+icomp);
796  Real delta = fine_arr(i ,j ,k,n) - d;
797  Real delta_xp = fine_arr(i+1,j ,k,n) - d_ip1;
798  Real delta_xm = fine_arr(i-1,j ,k,n) - d_im1;
799  Real delta_yp = fine_arr(i ,j+1,k,n) - d_jp1;
800  Real delta_ym = fine_arr(i ,j-1,k,n) - d_jm1;
801  Real Laplacian = delta_xp + delta_xm + delta_yp + delta_ym - 4.0*delta;
802  rhs_arr(i,j,k,n) += (F1*delta - F2*Laplacian) * Factor;
803  }
804  });
805  } // mfi
806  } // ivar_idx
807 }
void_bc
PhysBCFunctNoOp void_bc
Definition: ERF_InteriorGhostCells.cpp:5
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:143
IntVars::ymom
@ ymom
Definition: ERF_IndexDefines.H:141
IntVars::cons
@ cons
Definition: ERF_IndexDefines.H:139
IntVars::zmom
@ zmom
Definition: ERF_IndexDefines.H:142
IntVars::xmom
@ xmom
Definition: ERF_IndexDefines.H:140
Here is the call 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
121 {
122  BL_PROFILE_REGION("realbdy_compute_interior_ghost_RHS()");
123 
124  // NOTE: We pass the full width into this routine.
125  // For relaxation, the last cell is a halo
126  // cell for the Laplacian. We remove that
127  // cell here if it is present.
128 
129  // The width to do RHS augmentation
130  if (width > set_width+1) width -= 1;
131 
132  // Relaxation constants
133  Real F1 = 1./(10.*delta_t);
134  Real F2 = 1./(50.*delta_t);
135 
136  // Time interpolation
137  Real dT = bdy_time_interval;
138  Real time_since_start = time - start_bdy_time;
139  int n_time = static_cast<int>( time_since_start / dT);
140  Real alpha = (time_since_start - n_time * dT) / dT;
141  AMREX_ALWAYS_ASSERT( alpha >= 0. && alpha <= 1.0);
142  Real oma = 1.0 - alpha;
143 
144  /*
145  // UNIT TEST DEBUG
146  oma = 1.0; alpha = 0.0;
147  */
148 
149  // Temporary FABs for storage (owned/filled on all ranks)
150  FArrayBox U_xlo, U_xhi, U_ylo, U_yhi;
151  FArrayBox V_xlo, V_xhi, V_ylo, V_yhi;
152  FArrayBox T_xlo, T_xhi, T_ylo, T_yhi;
153 
154  // Variable index map (WRFBdyVars -> Vars)
155  Vector<int> var_map = {Vars::xvel, Vars::yvel, Vars::cons, Vars::cons};
156  Vector<int> ivar_map = {IntVars::xmom, IntVars::ymom, IntVars::cons, IntVars::cons};
157 
158  // Variable icomp map
159  Vector<int> comp_map = {0, 0, RhoTheta_comp};
160 
161  // Indices
162  int ivarU = RealBdyVars::U;
163  int ivarV = RealBdyVars::V;
164  int ivarT = RealBdyVars::T;
165  int BdyEnd = RealBdyVars::NumTypes-1;
166 
167 
168  // NOTE: These sizing of the temporary BDY FABS is
169  // GLOBAL and occurs over the entire BDY region.
170 
171  // Size the FABs
172  //==========================================================
173  for (int ivar(ivarU); ivar < BdyEnd; ivar++) {
174  int ivar_idx = var_map[ivar];
175  Box domain = geom.Domain();
176  auto ixtype = S_cur_data[ivar_idx].boxArray().ixType();
177  domain.convert(ixtype);
178 
179  // Grown domain to get the 4 halo boxes w/ ghost cells
180  // NOTE: 2 ghost cells needed here for Laplacian
181  // halo cell.
182  IntVect ng_vect{2,2,0};
183  Box gdom(domain); gdom.grow(ng_vect);
184  Box bx_xlo, bx_xhi, bx_ylo, bx_yhi;
185  compute_interior_ghost_bxs_xy(gdom, domain, width, 0,
186  bx_xlo, bx_xhi,
187  bx_ylo, bx_yhi,
188  ng_vect, true);
189 
190  // Size the FABs
191  if (ivar == ivarU) {
192  U_xlo.resize(bx_xlo,1,The_Async_Arena()); U_xhi.resize(bx_xhi,1,The_Async_Arena());
193  U_ylo.resize(bx_ylo,1,The_Async_Arena()); U_yhi.resize(bx_yhi,1,The_Async_Arena());
194  } else if (ivar == ivarV) {
195  V_xlo.resize(bx_xlo,1,The_Async_Arena()); V_xhi.resize(bx_xhi,1,The_Async_Arena());
196  V_ylo.resize(bx_ylo,1,The_Async_Arena()); V_yhi.resize(bx_yhi,1,The_Async_Arena());
197  } else if (ivar == ivarT){
198  T_xlo.resize(bx_xlo,1,The_Async_Arena()); T_xhi.resize(bx_xhi,1,The_Async_Arena());
199  T_ylo.resize(bx_ylo,1,The_Async_Arena()); T_yhi.resize(bx_yhi,1,The_Async_Arena());
200  } else {
201  continue;
202  }
203  } // ivar
204 
205 
206  // NOTE: These operations use the BDY FABS and RHO. The
207  // use of RHO to go from PRIM -> CONS requires that
208  // these operations be LOCAL. So we have allocated
209  // enough space to do global operations (1 rank) but
210  // will fill a subset of that data that the rank owns.
211 
212  // Populate FABs from bdy interpolation (primitive vars)
213  //==========================================================
214  for (int ivar(ivarU); ivar < BdyEnd; ivar++) {
215  int ivar_idx = var_map[ivar];
216  Box domain = geom.Domain();
217  auto ixtype = S_cur_data[ivar_idx].boxArray().ixType();
218  domain.convert(ixtype);
219  const auto& dom_lo = lbound(domain);
220  const auto& dom_hi = ubound(domain);
221 
222 #ifdef _OPENMP
223 #pragma omp parallel if (Gpu::notInLaunchRegion())
224 #endif
225  for ( MFIter mfi(S_cur_data[ivar_idx],TilingIfNotGPU()); mfi.isValid(); ++mfi) {
226  // NOTE: 2 ghost cells needed here for Laplacian
227  // halo cell.
228  IntVect ng_vect{2,2,0};
229  Box tbx = mfi.tilebox(ixtype.toIntVect(),ng_vect);
230  Box tbx_xlo, tbx_xhi, tbx_ylo, tbx_yhi;
231  compute_interior_ghost_bxs_xy(tbx, domain, width, 0,
232  tbx_xlo, tbx_xhi,
233  tbx_ylo, tbx_yhi,
234  ng_vect, true);
235 
236  Array4<Real> arr_xlo; Array4<Real> arr_xhi;
237  Array4<Real> arr_ylo; Array4<Real> arr_yhi;
238  if (ivar == ivarU) {
239  arr_xlo = U_xlo.array(); arr_xhi = U_xhi.array();
240  arr_ylo = U_ylo.array(); arr_yhi = U_yhi.array();
241  } else if (ivar == ivarV) {
242  arr_xlo = V_xlo.array(); arr_xhi = V_xhi.array();
243  arr_ylo = V_ylo.array(); arr_yhi = V_yhi.array();
244  } else if (ivar == ivarT){
245  arr_xlo = T_xlo.array(); arr_xhi = T_xhi.array();
246  arr_ylo = T_ylo.array(); arr_yhi = T_yhi.array();
247  } else {
248  continue;
249  }
250 
251  // Boundary data at fixed time intervals
252  const auto& bdatxlo_n = bdy_data_xlo[n_time ][ivar].const_array();
253  const auto& bdatxlo_np1 = bdy_data_xlo[n_time+1][ivar].const_array();
254  const auto& bdatxhi_n = bdy_data_xhi[n_time ][ivar].const_array();
255  const auto& bdatxhi_np1 = bdy_data_xhi[n_time+1][ivar].const_array();
256  const auto& bdatylo_n = bdy_data_ylo[n_time ][ivar].const_array();
257  const auto& bdatylo_np1 = bdy_data_ylo[n_time+1][ivar].const_array();
258  const auto& bdatyhi_n = bdy_data_yhi[n_time ][ivar].const_array();
259  const auto& bdatyhi_np1 = bdy_data_yhi[n_time+1][ivar].const_array();
260 
261  // Current density to convert to conserved vars
262  Array4<Real> r_arr = S_cur_data[IntVars::cons].array(mfi);
263 
264  // NOTE: width is now one less than the total bndy width
265  // if we have a relaxation zone; so we can access
266  // dom_lo/hi +- width. If we do not have a relax
267  // zone, this offset is set_width - 1.
268  int offset = set_width - 1;
269  if (width > set_width) offset = width;
270 
271  // Populate with interpolation (protect from ghost cells)
272  ParallelFor(tbx_xlo, tbx_xhi,
273  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
274  {
275  int ii = std::max(i , dom_lo.x);
276  ii = std::min(ii, dom_lo.x+offset);
277  int jj = std::max(j , dom_lo.y);
278  jj = std::min(jj, dom_hi.y);
279 
280  Real rho_interp;
281  if (ivar==ivarU) {
282  rho_interp = 0.5 * ( r_arr(i-1,j ,k) + r_arr(i,j,k) );
283  } else if (ivar==ivarV) {
284  rho_interp = 0.5 * ( r_arr(i ,j-1,k) + r_arr(i,j,k) );
285  } else {
286  rho_interp = r_arr(i,j,k);
287  }
288  arr_xlo(i,j,k) = rho_interp * ( oma * bdatxlo_n (ii,jj,k,0)
289  + alpha * bdatxlo_np1(ii,jj,k,0) );
290  },
291  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
292  {
293  int ii = std::max(i , dom_hi.x-offset);
294  ii = std::min(ii, dom_hi.x);
295  int jj = std::max(j , dom_lo.y);
296  jj = std::min(jj, dom_hi.y);
297 
298  Real rho_interp;
299  if (ivar==ivarU) {
300  rho_interp = 0.5 * ( r_arr(i-1,j ,k) + r_arr(i,j,k) );
301  } else if (ivar==ivarV) {
302  rho_interp = 0.5 * ( r_arr(i ,j-1,k) + r_arr(i,j,k) );
303  } else {
304  rho_interp = r_arr(i,j,k);
305  }
306  arr_xhi(i,j,k) = rho_interp * ( oma * bdatxhi_n (ii,jj,k,0)
307  + alpha * bdatxhi_np1(ii,jj,k,0) );
308  });
309 
310  ParallelFor(tbx_ylo, tbx_yhi,
311  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
312  {
313  int ii = std::max(i , dom_lo.x);
314  ii = std::min(ii, dom_hi.x);
315  int jj = std::max(j , dom_lo.y);
316  jj = std::min(jj, dom_lo.y+offset);
317 
318  Real rho_interp;
319  if (ivar==ivarU) {
320  rho_interp = 0.5 * ( r_arr(i-1,j ,k) + r_arr(i,j,k) );
321  } else if (ivar==ivarV) {
322  rho_interp = 0.5 * ( r_arr(i ,j-1,k) + r_arr(i,j,k) );
323  } else {
324  rho_interp = r_arr(i,j,k);
325  }
326  arr_ylo(i,j,k) = rho_interp * ( oma * bdatylo_n (ii,jj,k,0)
327  + alpha * bdatylo_np1(ii,jj,k,0) );
328  },
329  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
330  {
331  int ii = std::max(i , dom_lo.x);
332  ii = std::min(ii, dom_hi.x);
333  int jj = std::max(j , dom_hi.y-offset);
334  jj = std::min(jj, dom_hi.y);
335 
336  Real rho_interp;
337  if (ivar==ivarU) {
338  rho_interp = 0.5 * ( r_arr(i-1,j ,k) + r_arr(i,j,k) );
339  } else if (ivar==ivarV) {
340  rho_interp = 0.5 * ( r_arr(i ,j-1,k) + r_arr(i,j,k) );
341  } else {
342  rho_interp = r_arr(i,j,k);
343  }
344  arr_yhi(i,j,k) = rho_interp * ( oma * bdatyhi_n (ii,jj,k,0)
345  + alpha * bdatyhi_np1(ii,jj,k,0) );
346  });
347  } // mfi
348  } // ivar
349 
350 
351  // NOTE: These operations use current RHS, so they are
352  // LOCAL and occur over the data owned by a given rank.
353 
354  // Compute RHS in specified region
355  //==========================================================
356  if (set_width > 0) {
357  for (int ivar(ivarU); ivar < BdyEnd; ivar++) {
358  int ivar_idx = ivar_map[ivar];
359  int icomp = comp_map[ivar];
360 
361  Box domain = geom.Domain();
362  domain.convert(S_old_data[ivar_idx].boxArray().ixType());
363  const auto& dom_hi = ubound(domain);
364  const auto& dom_lo = lbound(domain);
365 
366 #ifdef _OPENMP
367 #pragma omp parallel if (Gpu::notInLaunchRegion())
368 #endif
369  for ( MFIter mfi(S_old_data[ivar_idx],TilingIfNotGPU()); mfi.isValid(); ++mfi) {
370  Box tbx = mfi.tilebox();
371  Box tbx_xlo, tbx_xhi, tbx_ylo, tbx_yhi;
372  compute_interior_ghost_bxs_xy(tbx, domain, width, 0,
373  tbx_xlo, tbx_xhi,
374  tbx_ylo, tbx_yhi);
375 
376  Array4<Real> rhs_arr; Array4<Real> data_arr;
377  Array4<Real> arr_xlo; Array4<Real> arr_xhi;
378  Array4<Real> arr_ylo; Array4<Real> arr_yhi;
379  if (ivar == ivarU) {
380  arr_xlo = U_xlo.array(); arr_xhi = U_xhi.array();
381  arr_ylo = U_ylo.array(); arr_yhi = U_yhi.array();
382  rhs_arr = S_rhs[IntVars::xmom].array(mfi);
383  data_arr = S_old_data[IntVars::xmom].array(mfi);
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  rhs_arr = S_rhs[IntVars::ymom].array(mfi);
388  data_arr = S_old_data[IntVars::ymom].array(mfi);
389  } else if (ivar == ivarT){
390  arr_xlo = T_xlo.array(); arr_xhi = T_xhi.array();
391  arr_ylo = T_ylo.array(); arr_yhi = T_yhi.array();
392  rhs_arr = S_rhs[IntVars::cons].array(mfi);
393  data_arr = S_old_data[IntVars::cons].array(mfi);
394  } else {
395  continue;
396  }
397 
398  realbdy_set_rhs_in_spec_region(delta_t, icomp, 1,
399  width, set_width, dom_lo, dom_hi,
400  tbx_xlo, tbx_xhi, tbx_ylo, tbx_yhi,
401  arr_xlo, arr_xhi, arr_ylo, arr_yhi,
402  data_arr, rhs_arr);
403  } // mfi
404  } // ivar
405  } // set_width
406 
407 
408  // NOTE: These operations use current density, so they are
409  // LOCAL and occur over the data owned by a given rank
410 
411  // Compute RHS in relaxation region
412  //==========================================================
413  if (width > set_width) {
414  for (int ivar(ivarU); ivar < BdyEnd; ivar++) {
415  int ivar_idx = ivar_map[ivar];
416  int icomp = comp_map[ivar];
417 
418  Box domain = geom.Domain();
419  domain.convert(S_cur_data[ivar_idx].boxArray().ixType());
420  const auto& dom_hi = ubound(domain);
421  const auto& dom_lo = lbound(domain);
422 
423  int width2 = width;
424  if (ivar_idx == IntVars::cons) width2 -= 1;
425 
426 #ifdef _OPENMP
427 #pragma omp parallel if (Gpu::notInLaunchRegion())
428 #endif
429  for ( MFIter mfi(S_cur_data[ivar_idx],TilingIfNotGPU()); mfi.isValid(); ++mfi) {
430  Box tbx = mfi.tilebox();
431  Box tbx_xlo, tbx_xhi, tbx_ylo, tbx_yhi;
432  compute_interior_ghost_bxs_xy(tbx, domain, width2, set_width,
433  tbx_xlo, tbx_xhi,
434  tbx_ylo, tbx_yhi);
435 
436  Array4<Real> rhs_arr; Array4<Real> data_arr;
437  Array4<Real> arr_xlo; Array4<Real> arr_xhi;
438  Array4<Real> arr_ylo; Array4<Real> arr_yhi;
439  if (ivar == ivarU) {
440  arr_xlo = U_xlo.array(); arr_xhi = U_xhi.array();
441  arr_ylo = U_ylo.array(); arr_yhi = U_yhi.array();
442  rhs_arr = S_rhs[IntVars::xmom].array(mfi);
443  data_arr = S_cur_data[IntVars::xmom].array(mfi);
444  } else if (ivar == ivarV) {
445  arr_xlo = V_xlo.array(); arr_xhi = V_xhi.array();
446  arr_ylo = V_ylo.array(); arr_yhi = V_yhi.array();
447  rhs_arr = S_rhs[IntVars::ymom].array(mfi);
448  data_arr = S_cur_data[IntVars::ymom].array(mfi);
449  } else if (ivar == ivarT){
450  arr_xlo = T_xlo.array(); arr_xhi = T_xhi.array();
451  arr_ylo = T_ylo.array(); arr_yhi = T_yhi.array();
452  rhs_arr = S_rhs[IntVars::cons].array(mfi);
453  data_arr = S_cur_data[IntVars::cons].array(mfi);
454  } else {
455  continue;
456  }
457 
458  realbdy_compute_laplacian_relaxation(icomp, 1,
459  width2, set_width, dom_lo, dom_hi, F1, F2,
460  tbx_xlo, tbx_xhi, tbx_ylo, tbx_yhi,
461  arr_xlo, arr_xhi, arr_ylo, arr_yhi,
462  data_arr, rhs_arr);
463 
464  /*
465  // UNIT TEST DEBUG
466  compute_interior_ghost_bxs_xy(tbx, domain, width+1, 0,
467  tbx_xlo, tbx_xhi,
468  tbx_ylo, tbx_yhi);
469  ParallelFor(tbx_xlo, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
470  {
471  if (arr_xlo(i,j,k) != data_arr(i,j,k,icomp)) {
472  Print() << "ERROR XLO: " << ivar << ' ' << icomp << ' ' << IntVect(i,j,k) << "\n";
473  exit(0);
474  }
475  });
476  ParallelFor(tbx_xhi, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
477  {
478  if (arr_xhi(i,j,k) != data_arr(i,j,k,icomp)) {
479  Print() << "ERROR XHI: " << ivar << ' ' << icomp << ' ' << IntVect(i,j,k) << "\n";
480  exit(0);
481  }
482  });
483  ParallelFor(tbx_ylo, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
484  {
485  if (arr_ylo(i,j,k) != data_arr(i,j,k,icomp)) {
486  Print() << "ERROR YLO: " << ivar << ' ' << icomp << ' ' << IntVect(i,j,k) << "\n";
487  exit(0);
488  }
489  });
490  ParallelFor(tbx_yhi, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
491  {
492  if (arr_yhi(i,j,k) != data_arr(i,j,k,icomp)) {
493  Print() << "ERROR YHI: " << ivar << ' ' << icomp << ' ' << IntVect(i,j,k) << "\n";
494  exit(0);
495  }
496  });
497  */
498  } // mfi
499  } // ivar
500  } // width
501 }
RhoTheta_comp
#define RhoTheta_comp
Definition: ERF_IndexDefines.H:37
compute_interior_ghost_bxs_xy
void compute_interior_ghost_bxs_xy(const Box &bx, const Box &domain, const int &width, const int &set_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_set_rhs_in_spec_region
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, 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:148
realbdy_compute_laplacian_relaxation
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:238
RealBdyVars::U
@ U
Definition: ERF_IndexDefines.H:97
RealBdyVars::NumTypes
@ NumTypes
Definition: ERF_IndexDefines.H:101
RealBdyVars::T
@ T
Definition: ERF_IndexDefines.H:99
RealBdyVars::V
@ V
Definition: ERF_IndexDefines.H:98
Vars::xvel
@ xvel
Definition: ERF_IndexDefines.H:130
Vars::cons
@ cons
Definition: ERF_IndexDefines.H:129
Vars::yvel
@ yvel
Definition: ERF_IndexDefines.H:131
Here is the call graph for this function:

Variable Documentation

◆ void_bc

PhysBCFunctNoOp void_bc

Referenced by fine_compute_interior_ghost_rhs().