ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
ERFPhysBCFunct_v Class Reference

#include <ERF_PhysBCFunct.H>

Collaboration diagram for ERFPhysBCFunct_v:

Public Member Functions

 ERFPhysBCFunct_v (const int lev, const amrex::Geometry &geom, const amrex::Vector< amrex::BCRec > &domain_bcs_type, const amrex::Gpu::DeviceVector< amrex::BCRec > &domain_bcs_type_d, amrex::Array< amrex::Array< amrex::Real, AMREX_SPACEDIM *2 >, AMREX_SPACEDIM+NBCVAR_max > bc_extdir_vals, amrex::Array< amrex::Array< amrex::Real, AMREX_SPACEDIM *2 >, AMREX_SPACEDIM+NBCVAR_max > bc_neumann_vals, std::unique_ptr< amrex::MultiFab > &z_phys_nd, const bool use_real_bcs, amrex::Real *v_bc_data)
 
 ~ERFPhysBCFunct_v ()
 
void operator() (amrex::MultiFab &mf, amrex::MultiFab &xvel, amrex::MultiFab &yvel, amrex::IntVect const &nghost, const double time, int bccomp, bool do_fb)
 
void impose_lateral_yvel_bcs (const amrex::Array4< amrex::Real > &dest_arr, const amrex::Array4< amrex::Real const > &xvel_arr, const amrex::Array4< amrex::Real const > &yvel_arr, const amrex::Box &bx, const amrex::Box &domain, int bccomp, const double time)
 
void impose_vertical_yvel_bcs (const amrex::Array4< amrex::Real > &dest_arr, const amrex::Box &bx, const amrex::Box &domain, const amrex::Array4< amrex::Real const > &z_nd, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > dxInv, int bccomp, const double time)
 

Private Attributes

int m_lev
 
amrex::Geometry m_geom
 
amrex::Vector< amrex::BCRec > m_domain_bcs_type
 
amrex::Gpu::DeviceVector< amrex::BCRec > m_domain_bcs_type_d
 
amrex::Array< amrex::Array< amrex::Real, AMREX_SPACEDIM *2 >, AMREX_SPACEDIM+NBCVAR_maxm_bc_extdir_vals
 
amrex::Array< amrex::Array< amrex::Real, AMREX_SPACEDIM *2 >, AMREX_SPACEDIM+NBCVAR_maxm_bc_neumann_vals
 
amrex::MultiFab * m_z_phys_nd
 
bool m_use_real_bcs
 
amrex::Realm_v_bc_data
 

Constructor & Destructor Documentation

◆ ERFPhysBCFunct_v()

ERFPhysBCFunct_v::ERFPhysBCFunct_v ( const int  lev,
const amrex::Geometry &  geom,
const amrex::Vector< amrex::BCRec > &  domain_bcs_type,
const amrex::Gpu::DeviceVector< amrex::BCRec > &  domain_bcs_type_d,
amrex::Array< amrex::Array< amrex::Real, AMREX_SPACEDIM *2 >, AMREX_SPACEDIM+NBCVAR_max bc_extdir_vals,
amrex::Array< amrex::Array< amrex::Real, AMREX_SPACEDIM *2 >, AMREX_SPACEDIM+NBCVAR_max bc_neumann_vals,
std::unique_ptr< amrex::MultiFab > &  z_phys_nd,
const bool  use_real_bcs,
amrex::Real v_bc_data 
)
inline

Construct the y-velocity physical boundary-condition functor.

Parameters
[in]levlevel index
[in]geomgeometry for the level
[in]domain_bcs_typehost boundary-condition records
[in]domain_bcs_type_ddevice boundary-condition records
[in]bc_extdir_valsexternal Dirichlet values for boundary faces
[in]bc_neumann_valsNeumann values for boundary faces
[in]z_phys_ndnodal physical height MultiFab
[in]use_real_bcswhether real-data boundary conditions are active
[in]v_bc_dataoptional y-velocity boundary data
257  : m_lev(lev),
258  m_geom(geom), m_domain_bcs_type(domain_bcs_type),
259  m_domain_bcs_type_d(domain_bcs_type_d),
260  m_bc_extdir_vals(bc_extdir_vals),
261  m_bc_neumann_vals(bc_neumann_vals),
262  m_z_phys_nd(z_phys_nd.get()),
263  m_use_real_bcs(use_real_bcs),
264  m_v_bc_data(v_bc_data)
265  { }
amrex::Geometry m_geom
Definition: ERF_PhysBCFunct.H:324
amrex::Gpu::DeviceVector< amrex::BCRec > m_domain_bcs_type_d
Definition: ERF_PhysBCFunct.H:326
amrex::Real * m_v_bc_data
Definition: ERF_PhysBCFunct.H:331
int m_lev
Definition: ERF_PhysBCFunct.H:323
amrex::Array< amrex::Array< amrex::Real, AMREX_SPACEDIM *2 >, AMREX_SPACEDIM+NBCVAR_max > m_bc_neumann_vals
Definition: ERF_PhysBCFunct.H:328
bool m_use_real_bcs
Definition: ERF_PhysBCFunct.H:330
amrex::Vector< amrex::BCRec > m_domain_bcs_type
Definition: ERF_PhysBCFunct.H:325
amrex::Array< amrex::Array< amrex::Real, AMREX_SPACEDIM *2 >, AMREX_SPACEDIM+NBCVAR_max > m_bc_extdir_vals
Definition: ERF_PhysBCFunct.H:327
amrex::MultiFab * m_z_phys_nd
Definition: ERF_PhysBCFunct.H:329

◆ ~ERFPhysBCFunct_v()

ERFPhysBCFunct_v::~ERFPhysBCFunct_v ( )
inline

Destroy the y-velocity physical boundary-condition functor.

270 {}

Member Function Documentation

◆ impose_lateral_yvel_bcs()

void ERFPhysBCFunct_v::impose_lateral_yvel_bcs ( const amrex::Array4< amrex::Real > &  dest_arr,
const amrex::Array4< amrex::Real const > &  xvel_arr,
const amrex::Array4< amrex::Real const > &  yvel_arr,
const amrex::Box &  bx,
const amrex::Box &  domain,
int  bccomp,
const double  time 
)

Impose lateral boundary conditions on y velocity.

Parameters
[in,out]dest_arry-velocity data to be filled
[in]xvel_arrx-velocity used to determine upwind inflow
[in]yvel_arry-velocity used to determine upwind inflow
[in]bxbox holding data to be filled
[in]domainsimulation domain
[in]bccompboundary-condition component
[in]timetime at which data should be filled

Impose lateral boundary conditions on y-component of velocity

Parameters
[in,out]dest_arrArray4 of the quantity to be filled
[in]xvel_arrx-velocity used to determine upwind Dirichlet inflow
[in]yvel_arry-velocity used to determine upwind Dirichlet inflow
[in]bxbox associated with this data
[in]domaincomputational domain
[in]bccompindex into m_domain_bcs_type
[in]timetime at which the data should be filled
22 {
23  BL_PROFILE_VAR("impose_lateral_yvel_bcs()",impose_lateral_yvel_bcs);
24  const auto& dom_lo = lbound(domain);
25  const auto& dom_hi = ubound(domain);
26 
27  // Based on BCRec for the domain, we need to make BCRec for this Box
28  // bccomp is used as starting index for m_domain_bcs_type
29  // 0 is used as starting index for bcrs
30  Vector<BCRec> bcrs(1);
31  setBC(enclosedCells(bx), domain, bccomp, 0, 1, m_domain_bcs_type, bcrs);
32 
33  // xlo: ori = 0
34  // ylo: ori = 1
35  // zlo: ori = 2
36  // xhi: ori = 3
37  // yhi: ori = 4
38  // zhi: ori = 5
39 
40  Gpu::DeviceVector<BCRec> bcrs_d(1);
41  Gpu::copyAsync(Gpu::hostToDevice, bcrs.begin(), bcrs.end(), bcrs_d.begin());
42  const BCRec* bc_ptr = bcrs_d.data();
43 
44  GpuArray<GpuArray<Real, AMREX_SPACEDIM*2>, 1> l_bc_extdir_vals_d;
45 
46  for (int ori = 0; ori < 2*AMREX_SPACEDIM; ori++) {
47  l_bc_extdir_vals_d[0][ori] = m_bc_extdir_vals[bccomp][ori];
48  }
49 
50  GeometryData const& geomdata = m_geom.data();
51  bool is_periodic_in_x = geomdata.isPeriodic(0);
52  bool is_periodic_in_y = geomdata.isPeriodic(1);
53 
54  // First do all ext_dir bcs
55  if (!is_periodic_in_x)
56  {
57  // Populate ghost cells on lo-x and hi-x domain boundaries
58  Real* yvel_bc_ptr = m_v_bc_data;
59  Box bx_xlo(bx); bx_xlo.setBig (0,dom_lo.x-1);
60  Box bx_xhi(bx); bx_xhi.setSmall(0,dom_hi.x+1);
61  ParallelFor(bx_xlo, bx_xhi,
62  [=] AMREX_GPU_DEVICE (int i, int j, int k) {
63  int iflip = dom_lo.x - 1- i;
64  if (bc_ptr[0].lo(0) == ERFBCType::ext_dir) {
65  dest_arr(i,j,k) = (yvel_bc_ptr) ? yvel_bc_ptr[k] : l_bc_extdir_vals_d[0][0];
66  } else if (bc_ptr[0].lo(0) == ERFBCType::ext_dir_upwind && xvel_arr(dom_lo.x,j,k) >= zero) {
67  dest_arr(i,j,k) = (yvel_bc_ptr) ? yvel_bc_ptr[k] : l_bc_extdir_vals_d[0][0];
68  } else if (bc_ptr[0].lo(0) == ERFBCType::foextrap) {
69  dest_arr(i,j,k) = dest_arr(dom_lo.x,j,k);
70  } else if (bc_ptr[0].lo(0) == ERFBCType::open) {
71  dest_arr(i,j,k) = dest_arr(dom_lo.x,j,k);
72  } else if (bc_ptr[0].lo(0) == ERFBCType::reflect_even) {
73  dest_arr(i,j,k) = dest_arr(iflip,j,k);
74  } else if (bc_ptr[0].lo(0) == ERFBCType::reflect_odd) {
75  dest_arr(i,j,k) = -dest_arr(iflip,j,k);
76  } else if (bc_ptr[0].lo(0) == ERFBCType::hoextrap) {
77  Real delta_i = (dom_lo.x - i);
78  dest_arr(i,j,k) = (one + delta_i)*dest_arr(dom_lo.x,j,k) - delta_i*dest_arr(dom_lo.x+1,j,k);
79  }
80  },
81  [=] AMREX_GPU_DEVICE (int i, int j, int k) {
82  int iflip = 2*dom_hi.x + 1 - i;
83  if (bc_ptr[0].hi(0) == ERFBCType::ext_dir) {
84  dest_arr(i,j,k) = (yvel_bc_ptr) ? yvel_bc_ptr[k] : l_bc_extdir_vals_d[0][3];
85  } else if (bc_ptr[0].hi(0) == ERFBCType::ext_dir_upwind && xvel_arr(dom_hi.x+1,j,k) <= zero) {
86  dest_arr(i,j,k) = (yvel_bc_ptr) ? yvel_bc_ptr[k] : l_bc_extdir_vals_d[0][3];
87  } else if (bc_ptr[0].hi(0) == ERFBCType::foextrap) {
88  dest_arr(i,j,k) = dest_arr(dom_hi.x,j,k);
89  } else if (bc_ptr[0].hi(0) == ERFBCType::open) {
90  dest_arr(i,j,k) = dest_arr(dom_hi.x,j,k);
91  } else if (bc_ptr[0].hi(0) == ERFBCType::reflect_even) {
92  dest_arr(i,j,k) = dest_arr(iflip,j,k);
93  } else if (bc_ptr[0].hi(0) == ERFBCType::reflect_odd) {
94  dest_arr(i,j,k) = -dest_arr(iflip,j,k);
95  } else if (bc_ptr[0].hi(0) == ERFBCType::hoextrap) {
96  Real delta_i = (i - dom_hi.x);
97  dest_arr(i,j,k) = (one + delta_i)*dest_arr(dom_hi.x,j,k) - delta_i*dest_arr(dom_hi.x-1,j,k);
98  }
99  }
100  );
101  }
102 
103  if (!is_periodic_in_y)
104  {
105  // Populate ghost cells on lo-y and hi-y domain boundaries
106  Real* yvel_bc_ptr = m_v_bc_data;
107  Box bx_ylo(bx); bx_ylo.setBig (1,dom_lo.y-1);
108  Box bx_yhi(bx); bx_yhi.setSmall(1,dom_hi.y+2);
109  Box bx_ylo_face(bx); bx_ylo_face.setSmall(1,dom_lo.y ); bx_ylo_face.setBig(1,dom_lo.y );
110  Box bx_yhi_face(bx); bx_yhi_face.setSmall(1,dom_hi.y+1); bx_yhi_face.setBig(1,dom_hi.y+1);
111  ParallelFor(bx_ylo, bx_ylo_face,
112  [=] AMREX_GPU_DEVICE (int i, int j, int k)
113  {
114  int jflip = dom_lo.y-j;
115  if (bc_ptr[0].lo(1) == ERFBCType::ext_dir) {
116  dest_arr(i,j,k) = (yvel_bc_ptr) ? yvel_bc_ptr[k] : l_bc_extdir_vals_d[0][1];
117  } else if (bc_ptr[0].lo(1) == ERFBCType::ext_dir_upwind && yvel_arr(i,dom_lo.y,k) >= zero) {
118  dest_arr(i,j,k) = (yvel_bc_ptr) ? yvel_bc_ptr[k] : l_bc_extdir_vals_d[0][1];
119  } else if (bc_ptr[0].lo(1) == ERFBCType::foextrap) {
120  dest_arr(i,j,k) = dest_arr(i,dom_lo.y,k);
121  } else if (bc_ptr[0].lo(1) == ERFBCType::open) {
122  dest_arr(i,j,k) = dest_arr(i,dom_lo.y,k);
123  } else if (bc_ptr[0].lo(1) == ERFBCType::reflect_even) {
124  dest_arr(i,j,k) = dest_arr(i,jflip,k);
125  } else if (bc_ptr[0].lo(1) == ERFBCType::reflect_odd) {
126  dest_arr(i,j,k) = -dest_arr(i,jflip,k);
127  } else if (bc_ptr[0].lo(1) == ERFBCType::neumann_int) {
128  dest_arr(i,j,k) = (Real(4.0)*dest_arr(i,dom_lo.y+1,k) - dest_arr(i,dom_lo.y+2,k))/three;
129  } else if (bc_ptr[0].lo(1) == ERFBCType::hoextrap) {
130  Real delta_j = (dom_lo.y - j);
131  dest_arr(i,j,k) = (one + delta_j)*dest_arr(i,dom_lo.y,k) - delta_j*dest_arr(i,dom_lo.y+1,k);
132  }
133  },
134  // We only set the values on the domain faces themselves if EXT_DIR or neumann_int
135  [=] AMREX_GPU_DEVICE (int i, int j, int k)
136  {
137  if (bc_ptr[0].lo(1) == ERFBCType::ext_dir) {
138  dest_arr(i,j,k) = (yvel_bc_ptr) ? yvel_bc_ptr[k] : l_bc_extdir_vals_d[0][1];
139  } else if (bc_ptr[0].lo(1) == ERFBCType::ext_dir_upwind && yvel_arr(i,dom_lo.y,k) >= zero) {
140  dest_arr(i,j,k) = (yvel_bc_ptr) ? yvel_bc_ptr[k] : l_bc_extdir_vals_d[0][1];
141  } else if (bc_ptr[0].lo(1) == ERFBCType::neumann_int) {
142  dest_arr(i,j,k) = (Real(4.0)*dest_arr(i,dom_lo.y+1,k) - dest_arr(i,dom_lo.y+2,k))/three;
143  }
144  }
145  );
146  ParallelFor(bx_yhi, bx_yhi_face,
147  [=] AMREX_GPU_DEVICE (int i, int j, int k)
148  {
149  int jflip = 2*(dom_hi.y + 1) - j;
150  if (bc_ptr[0].hi(1) == ERFBCType::ext_dir) {
151  dest_arr(i,j,k) = (yvel_bc_ptr) ? yvel_bc_ptr[k] : l_bc_extdir_vals_d[0][4];
152  } else if (bc_ptr[0].hi(1) == ERFBCType::ext_dir_upwind && yvel_arr(i,dom_hi.y+1,k) <= zero) {
153  dest_arr(i,j,k) = (yvel_bc_ptr) ? yvel_bc_ptr[k] : l_bc_extdir_vals_d[0][4];
154  } else if (bc_ptr[0].hi(1) == ERFBCType::foextrap) {
155  dest_arr(i,j,k) = dest_arr(i,dom_hi.y+1,k);
156  } else if (bc_ptr[0].hi(1) == ERFBCType::open) {
157  dest_arr(i,j,k) = dest_arr(i,dom_hi.y+1,k);
158  } else if (bc_ptr[0].hi(1) == ERFBCType::reflect_even) {
159  dest_arr(i,j,k) = dest_arr(i,jflip,k);
160  } else if (bc_ptr[0].hi(1) == ERFBCType::reflect_odd) {
161  dest_arr(i,j,k) = -dest_arr(i,jflip,k);
162  } else if (bc_ptr[0].hi(1) == ERFBCType::neumann_int) {
163  dest_arr(i,j,k) = (Real(4.0)*dest_arr(i,dom_hi.y,k) - dest_arr(i,dom_hi.y-1,k))/three;
164  } else if (bc_ptr[0].hi(1) == ERFBCType::hoextrap) {
165  Real delta_j = (j - dom_hi.y - 1);
166  dest_arr(i,j,k) = (one + delta_j)*dest_arr(i,dom_hi.y+1,k) - delta_j*dest_arr(i,dom_hi.y,k);
167  }
168  },
169  // We only set the values on the domain faces themselves if EXT_DIR or neumann_int
170  [=] AMREX_GPU_DEVICE (int i, int j, int k)
171  {
172  if (bc_ptr[0].hi(1) == ERFBCType::ext_dir) {
173  dest_arr(i,j,k) = (yvel_bc_ptr) ? yvel_bc_ptr[k] : l_bc_extdir_vals_d[0][4];
174  } else if (bc_ptr[0].hi(1) == ERFBCType::ext_dir_upwind && yvel_arr(i,dom_hi.y+1,k) <= zero) {
175  dest_arr(i,j,k) = (yvel_bc_ptr) ? yvel_bc_ptr[k] : l_bc_extdir_vals_d[0][4];
176  } else if (bc_ptr[0].hi(1) == ERFBCType::neumann_int) {
177  dest_arr(i,j,k) = (Real(4.0)*dest_arr(i,dom_hi.y,k) - dest_arr(i,dom_hi.y-1,k))/three;
178  }
179  }
180  );
181  }
182  Gpu::streamSynchronize();
183 }
constexpr amrex::Real three
Definition: ERF_Constants.H:11
constexpr amrex::Real one
Definition: ERF_Constants.H:9
constexpr amrex::Real zero
Definition: ERF_Constants.H:8
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);})
amrex::Real Real
Definition: ERF_ShocInterface.H:19
void impose_lateral_yvel_bcs(const amrex::Array4< amrex::Real > &dest_arr, const amrex::Array4< amrex::Real const > &xvel_arr, const amrex::Array4< amrex::Real const > &yvel_arr, const amrex::Box &bx, const amrex::Box &domain, int bccomp, const double time)
Definition: ERF_BoundaryConditionsYvel.cpp:17
@ open
Definition: ERF_IndexDefines.H:256
@ reflect_odd
Definition: ERF_IndexDefines.H:245
@ hoextrap
Definition: ERF_IndexDefines.H:250
@ foextrap
Definition: ERF_IndexDefines.H:248
@ ext_dir
Definition: ERF_IndexDefines.H:249
@ ext_dir_upwind
Definition: ERF_IndexDefines.H:257
@ neumann_int
Definition: ERF_IndexDefines.H:255
@ reflect_even
Definition: ERF_IndexDefines.H:247
Here is the call graph for this function:

◆ impose_vertical_yvel_bcs()

void ERFPhysBCFunct_v::impose_vertical_yvel_bcs ( const amrex::Array4< amrex::Real > &  dest_arr,
const amrex::Box &  bx,
const amrex::Box &  domain,
const amrex::Array4< amrex::Real const > &  z_nd,
const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM >  dxInv,
int  bccomp,
const double  time 
)

Impose vertical boundary conditions on y velocity.

Parameters
[in,out]dest_arry-velocity data to be filled
[in]bxbox holding data to be filled
[in]domainsimulation domain
[in]z_ndnodal physical height
[in]dxInvinverse cell size
[in]bccompboundary-condition component
[in]timetime at which data should be filled

Impose vertical boundary conditions on y-component of velocity

Parameters
[in,out]dest_arrthe Array4 of the quantity to be filled
[in]bxthe box associated with this data
[in]domainthe computational domain
[in]z_phys_ndheight coordinate at nodes
[in]dxInvinverse cell size array
[in]bccompindex into m_domain_bcs_type
[in]timetime at which the data should be filled
202 {
203  BL_PROFILE_VAR("impose_vertical_yvel_bcs()",impose_vertical_yvel_bcs);
204  const auto& dom_lo = lbound(domain);
205  const auto& dom_hi = ubound(domain);
206 
207  Box per_grown_domain(domain);
208  int growx = (m_geom.isPeriodic(0)) ? 1 : 0;
209  int growy = (m_geom.isPeriodic(1)) ? 1 : 0;
210  per_grown_domain.grow(IntVect(growx,growy,0));
211  const auto& perdom_lo = lbound(per_grown_domain);
212  const auto& perdom_hi = ubound(per_grown_domain);
213 
214  // Based on BCRec for the domain, we need to make BCRec for this Box
215  // bccomp is used as starting index for m_domain_bcs_type
216  // 0 is used as starting index for bcrs
217  Vector<BCRec> bcrs(1);
218  setBC(enclosedCells(bx), domain, bccomp, 0, 1, m_domain_bcs_type, bcrs);
219 
220  // xlo: ori = 0
221  // ylo: ori = 1
222  // zlo: ori = 2
223  // xhi: ori = 3
224  // yhi: ori = 4
225  // zhi: ori = 5
226 
227  Gpu::DeviceVector<BCRec> bcrs_d(1);
228  Gpu::copyAsync(Gpu::hostToDevice, bcrs.begin(), bcrs.end(), bcrs_d.begin());
229  const BCRec* bc_ptr = bcrs_d.data();
230 
231  GpuArray<GpuArray<Real, AMREX_SPACEDIM*2>, 1> l_bc_extdir_vals_d;
232 
233  for (int ori = 0; ori < 2*AMREX_SPACEDIM; ori++) {
234  l_bc_extdir_vals_d[0][ori] = m_bc_extdir_vals[bccomp][ori];
235  }
236 
237  GeometryData const& geomdata = m_geom.data();
238 
239  {
240  // Populate ghost cells on lo-z and hi-z domain boundaries
241  Box bx_zlo(bx); bx_zlo.setBig (2,dom_lo.z-1);
242  Box bx_zhi(bx); bx_zhi.setSmall(2,dom_hi.z+1);
243  ParallelFor(bx_zlo, bx_zhi,
244  [=] AMREX_GPU_DEVICE (int i, int j, int k) {
245  int kflip = dom_lo.z - 1 - k;
246  if (bc_ptr[0].lo(2) == ERFBCType::ext_dir) {
247  dest_arr(i,j,k) = l_bc_extdir_vals_d[0][2];
248  } else if (bc_ptr[0].lo(2) == ERFBCType::foextrap) {
249  dest_arr(i,j,k) = dest_arr(i,j,dom_lo.z);
250  } else if (bc_ptr[0].lo(2) == ERFBCType::open) {
251  dest_arr(i,j,k) = dest_arr(i,j,dom_lo.z);
252  } else if (bc_ptr[0].lo(2) == ERFBCType::reflect_even) {
253  dest_arr(i,j,k) = dest_arr(i,j,kflip);
254  } else if (bc_ptr[0].lo(2) == ERFBCType::reflect_odd) {
255  dest_arr(i,j,k) = -dest_arr(i,j,kflip);
256  } else if (bc_ptr[0].lo(2) == ERFBCType::hoextrap) {
257  Real delta_k = (dom_lo.z - k);
258  dest_arr(i,j,k) = (one + delta_k)*dest_arr(i,j,dom_lo.z) - delta_k*dest_arr(i,j,dom_lo.z+1);
259  }
260  },
261  [=] AMREX_GPU_DEVICE (int i, int j, int k) {
262  int kflip = 2*dom_hi.z + 1 - k;
263  if (bc_ptr[0].hi(2) == ERFBCType::ext_dir) {
264  dest_arr(i,j,k) = l_bc_extdir_vals_d[0][5];
265  } else if (bc_ptr[0].hi(2) == ERFBCType::foextrap) {
266  dest_arr(i,j,k) = dest_arr(i,j,dom_hi.z);
267  } else if (bc_ptr[0].hi(2) == ERFBCType::open) {
268  dest_arr(i,j,k) = dest_arr(i,j,dom_hi.z);
269  } else if (bc_ptr[0].hi(2) == ERFBCType::reflect_even) {
270  dest_arr(i,j,k) = dest_arr(i,j,kflip);
271  } else if (bc_ptr[0].hi(2) == ERFBCType::reflect_odd) {
272  dest_arr(i,j,k) = -dest_arr(i,j,kflip);
273  } else if (bc_ptr[0].hi(2) == ERFBCType::hoextrap){
274  Real delta_k = (k - dom_hi.z);
275  dest_arr(i,j,k) = (one + delta_k)*dest_arr(i,j,dom_hi.z) - delta_k*dest_arr(i,j,dom_hi.z-1);
276  }
277  }
278  );
279  }
280 
281  if (m_z_phys_nd) {
282 
283  const auto& bx_lo = lbound(bx);
284  const auto& bx_hi = ubound(bx);
285 
286  const auto& zphys_lo = lbound(Box(z_phys_nd));
287  const auto& zphys_hi = ubound(Box(z_phys_nd));
288 
289  // Neumann conditions (d<var>/dn = 0) must be aware of the surface normal with terrain.
290  // An additional source term arises from d<var>/dx & d<var>/dy & met_h_xi/eta/zeta.
291  //=====================================================================================
292  // Only modify scalars, U, or V
293  // Loop over each component
294  // Hit for Neumann condition at kmin
295  if(bcrs[0].lo(2) == ERFBCType::foextrap) {
296  // Loop over ghost cells in bottom XY-plane (valid box)
297  Box xybx = bx;
298  xybx.setBig(2,-1);
299  xybx.setSmall(2,bx.smallEnd()[2]);
300  int k0 = 0;
301 
302  // Get the dz cell size
303  Real dz = geomdata.CellSize(2);
304 
305  // Fill all the Neumann srcs with terrain
306  ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
307  {
308  // Clip indices for ghost-cells
309  int ii = amrex::min(amrex::max(i,perdom_lo.x),perdom_hi.x);
310  ii = amrex::min(amrex::max(ii,zphys_lo.x),zphys_hi.x);
311  int jj = amrex::min(amrex::max(j,perdom_lo.y),perdom_hi.y);
312  jj = amrex::min(amrex::max(jj,zphys_lo.y),zphys_hi.y);
313 
314  // Get metrics
315  Real met_h_xi = Compute_h_xi_AtJface (ii, jj, k0, dxInv, z_phys_nd);
316  Real met_h_eta = Compute_h_eta_AtJface (ii, jj, k0, dxInv, z_phys_nd);
317  Real met_h_zeta = Compute_h_zeta_AtJface(ii, jj, k0, dxInv, z_phys_nd);
318 
319  // GradX at IJK location inside domain -- this relies on the assumption that we have
320  // used foextrap for cell-centered quantities outside the domain to define the gradient as zero
321  Real GradVarx, GradVary;
322  if ( i < dom_lo.x-1 || i > dom_hi.x+1 || (i+1 > bx_hi.x && i-1 < bx_lo.x) ) {
323  GradVarx = zero;
324  } else if (i+1 > bx_hi.x) {
325  GradVarx = dxInv[0] * (dest_arr(i ,j,k0) - dest_arr(i-1,j,k0));
326  } else if (i-1 < bx_lo.x) {
327  GradVarx = dxInv[0] * (dest_arr(i+1,j,k0) - dest_arr(i ,j,k0));
328  } else {
329  GradVarx = myhalf * dxInv[0] * (dest_arr(i+1,j,k0) - dest_arr(i-1,j,k0));
330  }
331 
332  // GradY at IJK location inside domain -- this relies on the assumption that we have
333  // used foextrap for cell-centered quantities outside the domain to define the gradient as zero
334  if ( j < dom_lo.y-1 || j > dom_hi.y+1 || (j+1 > bx_hi.y && j-1 < bx_lo.y) ) {
335  GradVary = zero;
336  } else if (j+1 > bx_hi.y) {
337  GradVary = dxInv[1] * (dest_arr(i,j ,k0) - dest_arr(i,j-1,k0));
338  } else if (j-1 < bx_lo.y) {
339  GradVary = dxInv[1] * (dest_arr(i,j+1,k0) - dest_arr(i,j ,k0));
340  } else {
341  GradVary = myhalf * dxInv[1] * (dest_arr(i,j+1,k0) - dest_arr(i,j-1,k0));
342  }
343 
344  // Prefactor
345  Real met_fac = met_h_zeta / ( met_h_xi*met_h_xi + met_h_eta*met_h_eta + one );
346 
347  // Accumulate in bottom ghost cell (EXTRAP already populated)
348  dest_arr(i,j,k) -= dz * met_fac * ( met_h_xi * GradVarx + met_h_eta * GradVary );
349  });
350  } // foextrap
351  } //m_z_phys_nd
352  Gpu::streamSynchronize();
353 }
constexpr amrex::Real myhalf
Definition: ERF_Constants.H:13
amrex::GpuArray< Real, AMREX_SPACEDIM > dxInv
Definition: ERF_InitCustomPertVels_ParticleTests.H:17
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_xi_AtJface(const int &i, const int &j, const int &k, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &z_nd)
Definition: ERF_TerrainMetrics.H:157
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_zeta_AtJface(const int &i, const int &j, const int &k, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &z_nd)
Definition: ERF_TerrainMetrics.H:144
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_eta_AtJface(const int &i, const int &j, const int &k, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &z_nd)
Definition: ERF_TerrainMetrics.H:170
void impose_vertical_yvel_bcs(const amrex::Array4< amrex::Real > &dest_arr, const amrex::Box &bx, const amrex::Box &domain, const amrex::Array4< amrex::Real const > &z_nd, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > dxInv, int bccomp, const double time)
Definition: ERF_BoundaryConditionsYvel.cpp:197
@ dz
Definition: ERF_AdvanceWSM6.cpp:104
Here is the call graph for this function:

◆ operator()()

void ERFPhysBCFunct_v::operator() ( amrex::MultiFab &  mf,
amrex::MultiFab &  xvel,
amrex::MultiFab &  yvel,
amrex::IntVect const &  nghost,
const double  time,
int  bccomp,
bool  do_fb 
)

Impose physical boundary conditions at domain boundaries

Parameters
[in,out]mfy-velocity MultiFab to be filled
[in]xvelx-velocity MultiFab used by upwind boundary conditions
[in]yvely-velocity MultiFab used by upwind boundary conditions
[in]nghostnumber of ghost cells to be filled for velocity components
[in]timetime at which the data should be filled
[in]bccompindex into m_domain_bcs_type
[in]do_fbwhether to fill periodic/interior ghost cells first

Impose physical boundary conditions on y velocity.

Parameters
[in,out]mfy-velocity MultiFab to be filled
[in]xvelx-velocity MultiFab used by upwind boundary conditions
[in]yvely-velocity MultiFab used by upwind boundary conditions
[in]nghostnumber of ghost cells to be filled
[in]timetime at which data should be filled
[in]bccompindex into m_domain_bcs_type
[in]do_fbwhether to fill periodic/interior ghost cells first
187 {
188  BL_PROFILE("ERFPhysBCFunct_v::()");
189 
190  //
191  // We fill all of the interior and periodic ghost cells first, so we can fill
192  // those directly inside the lateral and vertical calls.
193  // If triply periodic this is all we do
194  //
195  if (do_fb) {
196  mf.FillBoundary(m_geom.periodicity());
197  }
198 
199  if (m_geom.isAllPeriodic()) return;
200 
201  const auto& domain = m_geom.Domain();
202  const auto dxInv = m_geom.InvCellSizeArray();
203 
204  Box gdomainy = surroundingNodes(domain,1);
205  for (int i = 0; i < AMREX_SPACEDIM; ++i) {
206  if (m_geom.isPeriodic(i)) {
207  gdomainy.grow(i, nghost[i]);
208  }
209  }
210 
211 #ifdef AMREX_USE_OMP
212 #pragma omp parallel if (Gpu::notInLaunchRegion())
213 #endif
214  {
215  for (MFIter mfi(mf,false); mfi.isValid(); ++mfi)
216  {
217  //
218  // This is the box we pass to the different routines
219  // NOTE -- this is the full grid NOT the tile box
220  //
221  Box bx = mfi.validbox();
222 
223  //
224  // These are the boxes we use to test on relative to the domain
225  //
226  Box ybx1 = surroundingNodes(bx,1); ybx1.grow(nghost);
227  if (ybx1.smallEnd(2) < domain.smallEnd(2)) ybx1.setSmall(2,domain.smallEnd(2));
228  if (ybx1.bigEnd(2) > domain.bigEnd(2)) ybx1.setBig(2,domain.bigEnd(2));
229 
230  Box ybx2 = surroundingNodes(bx,1); ybx2.grow(nghost);
231 
232  Array4<const Real> z_nd_arr;
233 
234  if (m_z_phys_nd) {
235  z_nd_arr = m_z_phys_nd->const_array(mfi);
236  }
237 
238  if (!gdomainy.contains(ybx2)) {
239  Array4< Real> const& dest_arr = mf.array(mfi);
240  Array4<const Real> const& velx_arr = xvel.const_array(mfi);
241  Array4<const Real> const& vely_arr = yvel.const_array(mfi);
242 
243  if (!m_use_real_bcs) {
244  impose_lateral_yvel_bcs(dest_arr,velx_arr,vely_arr,ybx1,domain,bccomp,time);
245  }
246 
247  impose_vertical_yvel_bcs(dest_arr,ybx2,domain,z_nd_arr,dxInv,bccomp,time);
248  }
249 
250  } // MFIter
251  } // OpenMP
252 } // operator()
@ xvel
Definition: ERF_IndexDefines.H:177
@ yvel
Definition: ERF_IndexDefines.H:178

Member Data Documentation

◆ m_bc_extdir_vals

amrex::Array<amrex::Array<amrex::Real, AMREX_SPACEDIM*2>,AMREX_SPACEDIM+NBCVAR_max> ERFPhysBCFunct_v::m_bc_extdir_vals
private

◆ m_bc_neumann_vals

amrex::Array<amrex::Array<amrex::Real, AMREX_SPACEDIM*2>,AMREX_SPACEDIM+NBCVAR_max> ERFPhysBCFunct_v::m_bc_neumann_vals
private

◆ m_domain_bcs_type

amrex::Vector<amrex::BCRec> ERFPhysBCFunct_v::m_domain_bcs_type
private

◆ m_domain_bcs_type_d

amrex::Gpu::DeviceVector<amrex::BCRec> ERFPhysBCFunct_v::m_domain_bcs_type_d
private

◆ m_geom

amrex::Geometry ERFPhysBCFunct_v::m_geom
private

◆ m_lev

int ERFPhysBCFunct_v::m_lev
private

◆ m_use_real_bcs

bool ERFPhysBCFunct_v::m_use_real_bcs
private

◆ m_v_bc_data

amrex::Real* ERFPhysBCFunct_v::m_v_bc_data
private

◆ m_z_phys_nd

amrex::MultiFab* ERFPhysBCFunct_v::m_z_phys_nd
private

The documentation for this class was generated from the following files: