ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
Public Member Functions | Private Attributes | List of all members
ERFPhysBCFunct_cons Class Reference

#include <ERF_PhysBCFunct.H>

Collaboration diagram for ERFPhysBCFunct_cons:
Collaboration graph
[legend]

Public Member Functions

 ERFPhysBCFunct_cons (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, const TerrainType &terrain_type, std::unique_ptr< amrex::MultiFab > &z_phys_nd, const bool use_real_bcs, amrex::Real *th_bc_data, const amrex::Array< int, AMREX_SPACEDIM *2 > &th_file_face)
 
 ~ERFPhysBCFunct_cons ()
 
void operator() (amrex::MultiFab &mf, amrex::MultiFab &xvel, amrex::MultiFab &yvel, int icomp, int ncomp, amrex::IntVect const &nghost, const double time, int bccomp_cons, bool do_fb=true, bool do_terrain_adjustment=true)
 
void impose_lateral_cons_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 icomp, int ncomp, amrex::IntVect ng, const double time)
 
void impose_vertical_cons_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 icomp, int ncomp, const double time, bool do_terrain_adjustment=true)
 

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_max > m_bc_extdir_vals
 
amrex::Array< amrex::Array< amrex::Real, AMREX_SPACEDIM *2 >, AMREX_SPACEDIM+NBCVAR_max > m_bc_neumann_vals
 
TerrainType m_terrain_type
 
amrex::MultiFab * m_z_phys_nd
 
bool m_use_real_bcs
 
amrex::Real * m_th_bc_data
 
amrex::GpuArray< int, AMREX_SPACEDIM *2 > m_th_file_face
 
BelowGroundGhostSync m_below_ground_sync
 

Constructor & Destructor Documentation

◆ ERFPhysBCFunct_cons()

ERFPhysBCFunct_cons::ERFPhysBCFunct_cons ( 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,
const TerrainType &  terrain_type,
std::unique_ptr< amrex::MultiFab > &  z_phys_nd,
const bool  use_real_bcs,
amrex::Real *  th_bc_data,
const amrex::Array< int, AMREX_SPACEDIM *2 > &  th_file_face 
)
inline

Construct the conserved-variable 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]terrain_typeterrain type for the level
[in]z_phys_ndnodal physical height MultiFab
[in]use_real_bcswhether real-data boundary conditions are active
[in]th_bc_dataoptional theta boundary data (one profile in z)
[in]th_file_faceper face: 0 the face's theta is not th_bc_data, 1 th_bc_data holds rho*theta for that face, 2 it holds theta
83  : m_lev(lev), m_geom(geom),
84  m_domain_bcs_type(domain_bcs_type),
85  m_domain_bcs_type_d(domain_bcs_type_d),
86  m_bc_extdir_vals(bc_extdir_vals),
87  m_bc_neumann_vals(bc_neumann_vals),
88  m_terrain_type(terrain_type),
89  m_z_phys_nd(z_phys_nd.get()),
90  m_use_real_bcs(use_real_bcs),
91  m_th_bc_data(th_bc_data)
92  {
93  for (int ori = 0; ori < 2*AMREX_SPACEDIM; ++ori) { m_th_file_face[ori] = th_file_face[ori]; }
94  }
ERFPhysBCFunct_cons::m_use_real_bcs
bool m_use_real_bcs
Definition: ERF_PhysBCFunct.H:168
ERFPhysBCFunct_cons::m_geom
amrex::Geometry m_geom
Definition: ERF_PhysBCFunct.H:161
ERFPhysBCFunct_cons::m_bc_extdir_vals
amrex::Array< amrex::Array< amrex::Real, AMREX_SPACEDIM *2 >, AMREX_SPACEDIM+NBCVAR_max > m_bc_extdir_vals
Definition: ERF_PhysBCFunct.H:164
ERFPhysBCFunct_cons::m_bc_neumann_vals
amrex::Array< amrex::Array< amrex::Real, AMREX_SPACEDIM *2 >, AMREX_SPACEDIM+NBCVAR_max > m_bc_neumann_vals
Definition: ERF_PhysBCFunct.H:165
ERFPhysBCFunct_cons::m_domain_bcs_type_d
amrex::Gpu::DeviceVector< amrex::BCRec > m_domain_bcs_type_d
Definition: ERF_PhysBCFunct.H:163
ERFPhysBCFunct_cons::m_domain_bcs_type
amrex::Vector< amrex::BCRec > m_domain_bcs_type
Definition: ERF_PhysBCFunct.H:162
ERFPhysBCFunct_cons::m_terrain_type
TerrainType m_terrain_type
Definition: ERF_PhysBCFunct.H:166
ERFPhysBCFunct_cons::m_th_bc_data
amrex::Real * m_th_bc_data
Definition: ERF_PhysBCFunct.H:169
ERFPhysBCFunct_cons::m_z_phys_nd
amrex::MultiFab * m_z_phys_nd
Definition: ERF_PhysBCFunct.H:167
ERFPhysBCFunct_cons::m_lev
int m_lev
Definition: ERF_PhysBCFunct.H:160
ERFPhysBCFunct_cons::m_th_file_face
amrex::GpuArray< int, AMREX_SPACEDIM *2 > m_th_file_face
Definition: ERF_PhysBCFunct.H:170

◆ ~ERFPhysBCFunct_cons()

ERFPhysBCFunct_cons::~ERFPhysBCFunct_cons ( )
inline

Destroy the conserved-variable physical boundary-condition functor.

99 {}

Member Function Documentation

◆ impose_lateral_cons_bcs()

void ERFPhysBCFunct_cons::impose_lateral_cons_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  icomp,
int  ncomp,
amrex::IntVect  ng,
const double  time 
)

Impose lateral boundary conditions on conserved variables.

Parameters
[in,out]dest_arrcell-centered 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]icompstarting component
[in]ncompnumber of components
[in]ngnumber of ghost cells
[in]timetime at which data should be filled

Impose lateral boundary conditions on conserved scalars (at cell centers)

Parameters
[in,out]dest_arrcell-centered data 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 holding data to be filled
[in]domainsimulation domain
[in]icompindex into the MultiFab – this can be any value from 0 to NVAR-1
[in]ncompthe number of components – this can be any value from 1 to NVAR as long as icomp+ncomp <= NVAR-one
[in]ngnumber of ghost cells in each coordinate direction
[in]timetime at which the data should be filled
28 {
29  BL_PROFILE_VAR("impose_lateral_cons_bcs()",impose_lateral_cons_bcs);
30  const auto& dom_lo = lbound(domain);
31  const auto& dom_hi = ubound(domain);
32 
33  //
34  // An ext_dir_upwind face decides inflow against the face velocity, but the box we
35  // loop over does not always have a face velocity at every one of its cells: the
36  // state carries one more lateral ghost than the velocities do (ngrow_state =
37  // ComputeGhostCells()+1 against ngrow_vels = ComputeGhostCells()), and a velocity
38  // box is nodal in its own direction where the array being read is cell centered.
39  // Either way the outermost ghost cells have no velocity of their own, so sample
40  // the nearest one that exists. Cells already inside the velocity array are
41  // unaffected; only that outermost layer changes, and it previously read out of
42  // bounds.
43  //
44  const auto& xvel_lo = lbound(xvel_arr);
45  const auto& xvel_hi = ubound(xvel_arr);
46  const auto& yvel_lo = lbound(yvel_arr);
47  const auto& yvel_hi = ubound(yvel_arr);
48 
49  // xlo: ori = 0
50  // ylo: ori = 1
51  // zlo: ori = 2
52  // xhi: ori = 3
53  // yhi: ori = 4
54  // zhi: ori = 5
55 
56  // Based on BCRec for the domain, we need to make BCRec for this Box
57  // 0 is used as starting index for bcrs
58  Vector<BCRec> bcrs(ncomp);
59 
60  GpuArray<GpuArray<Real, AMREX_SPACEDIM*2>,NBCVAR_max> l_bc_extdir_vals_d;
61 
62  const int* bxlo = bx.loVect();
63  const int* bxhi = bx.hiVect();
64  const int* dlo = domain.loVect();
65  const int* dhi = domain.hiVect();
66 
67  for (int nc = 0; nc < ncomp; nc++)
68  {
69  int bc_comp = (icomp+nc >= RhoScalar_comp && icomp+nc < RhoScalar_comp+NSCALARS) ?
70  BCVars::RhoScalar_bc_comp : icomp+nc;
71  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
72  for (int dir = 0; dir < AMREX_SPACEDIM; dir++)
73  {
74  bcrs[nc].setLo(dir, ( bxlo[dir]<=dlo[dir]
75  ? m_domain_bcs_type[bc_comp].lo(dir) : BCType::int_dir ));
76  bcrs[nc].setHi(dir, ( bxhi[dir]>=dhi[dir]
77  ? m_domain_bcs_type[bc_comp].hi(dir) : BCType::int_dir ));
78  }
79 
80  for (int ori = 0; ori < 2*AMREX_SPACEDIM; ori++) {
81  l_bc_extdir_vals_d[bc_comp][ori] = m_bc_extdir_vals[bc_comp][ori];
82  }
83  }
84 
85  Gpu::DeviceVector<BCRec> bcrs_d(ncomp);
86  Gpu::copyAsync(Gpu::hostToDevice, bcrs.begin(), bcrs.end(), bcrs_d.begin());
87  const BCRec* bc_ptr = bcrs_d.data();
88 
89  GeometryData const& geomdata = m_geom.data();
90  bool is_periodic_in_x = geomdata.isPeriodic(0);
91  bool is_periodic_in_y = geomdata.isPeriodic(1);
92 
93  // The density a ghost cell on each face carries after this fill: the prescribed inflow
94  // density where Rho itself is ext_dir or ext_dir_upwind on that face (<face>.density was
95  // given), otherwise
96  // the first interior cell's. A primitive Dirichlet value (ext_dir_prim, e.g. theta from a
97  // dirichlet_file) must be multiplied by that density, not by the interior density alone:
98  // otherwise the ghost holds rho_interior * theta over a ghost density rho_in, and the
99  // primitive theta the scheme sees is scaled by rho_interior / rho_in.
100  // Only the faces that read a theta file use th_bc_data, in the form they read it in
101  // (one profile in z, so the lateral faces only); see ERF::m_th_file_face
102  const GpuArray<int,AMREX_SPACEDIM*2> th_file = m_th_file_face;
103  GpuArray<int,6> rho_is_extdir;
104  GpuArray<Real,6> rho_extdir_val;
105  for (int ori = 0; ori < 2*AMREX_SPACEDIM; ori++) {
106  const int dir = ori % AMREX_SPACEDIM;
107  const bool low = (ori < AMREX_SPACEDIM);
108  const int rho_type = low ? m_domain_bcs_type[BCVars::Rho_bc_comp].lo(dir)
109  : m_domain_bcs_type[BCVars::Rho_bc_comp].hi(dir);
110  const Real rho_val = m_bc_extdir_vals[BCVars::Rho_bc_comp][ori];
111  // ext_dir_upwind: the ghost takes rho_in where the face flows in, which is the only
112  // place a Dirichlet theta is applied on such a face
113  rho_is_extdir[ori] = ((rho_type == ERFBCType::ext_dir || rho_type == ERFBCType::ext_dir_upwind)
114  && rho_val > zero) ? 1 : 0;
115  rho_extdir_val[ori] = rho_val;
116  }
117  // Not covered: a density ingested from boundary planes (Rho typed ext_dir_ingested) varies
118  // over the face, so a primitive theta file on that face is still converted with the first
119  // interior cell's density below; the planes then carry the matching theta in practice.
120 
121  // First do all ext_dir bcs
122  if (!is_periodic_in_x)
123  {
124  Real* th_bc_ptr = m_th_bc_data;
125  Box bx_xlo(bx); bx_xlo.setBig (0,dom_lo.x-1);
126  Box bx_xhi(bx); bx_xhi.setSmall(0,dom_hi.x+1);
127  //
128  // If we are setting Dirichlet values, set them in all ghost cells on an inflow face
129  // "bx" is already grown in the x- and y-directions so here we just grow it in z
130  //
131  bx_xlo.grow(2,ng[2]);
132  bx_xhi.grow(2,ng[2]);
133  ParallelFor(
134  bx_xlo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n)
135  {
136  int dest_comp = icomp+n;
137  int k_profile = amrex::min(amrex::max(k, dom_lo.z), dom_hi.z) - dom_lo.z;
138  int bc_comp = (dest_comp >= RhoScalar_comp && dest_comp < RhoScalar_comp+NSCALARS) ?
139  BCVars::RhoScalar_bc_comp : dest_comp;
140  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
141  int l_bc_type = bc_ptr[n].lo(0);
142 
143  if ( (l_bc_type == ERFBCType::ext_dir) ||
144  (l_bc_type == ERFBCType::ext_dir_upwind && xvel_arr(dom_lo.x,Clamp(j,xvel_lo.y,xvel_hi.y),Clamp(k,xvel_lo.z,xvel_hi.z)) >= zero) )
145  {
146  if ((dest_comp == RhoTheta_comp) && th_bc_ptr && th_file[0] != 0) {
147  // a theta profile reaches ext_dir only as rho*theta (th_file 1) or on an
148  // ext_dir_upwind face, which has no primitive variant (th_file 2):
149  // convert the latter with the ghost density here
150  const Real rho_th = (th_file[0] == 2)
151  ? (rho_is_extdir[0] ? rho_extdir_val[0] : dest_arr(dom_lo.x,j,k,Rho_comp))
152  : Real(1.0);
153  dest_arr(i,j,k,dest_comp) = rho_th * th_bc_ptr[k_profile];
154  } else {
155  dest_arr(i,j,k,dest_comp) = l_bc_extdir_vals_d[bc_comp][0];
156  }
157  } else if (l_bc_type == ERFBCType::ext_dir_prim) {
158  Real rho = rho_is_extdir[0] ? rho_extdir_val[0] : dest_arr(dom_lo.x,j,k,Rho_comp);
159  if ((dest_comp == RhoTheta_comp) && th_bc_ptr && th_file[0] == 2) {
160  dest_arr(i,j,k,dest_comp) = rho * th_bc_ptr[k_profile];
161  } else {
162  dest_arr(i,j,k,dest_comp) = rho * l_bc_extdir_vals_d[bc_comp][0];
163  }
164  }
165  },
166  bx_xhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n)
167  {
168  int dest_comp = icomp+n;
169  int k_profile = amrex::min(amrex::max(k, dom_lo.z), dom_hi.z) - dom_lo.z;
170  int bc_comp = (dest_comp >= RhoScalar_comp && dest_comp < RhoScalar_comp+NSCALARS) ?
171  BCVars::RhoScalar_bc_comp : dest_comp;
172  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
173  int h_bc_type = bc_ptr[n].hi(0);
174 
175  if ( (h_bc_type == ERFBCType::ext_dir) ||
176  (h_bc_type == ERFBCType::ext_dir_upwind && xvel_arr(dom_hi.x+1,Clamp(j,xvel_lo.y,xvel_hi.y),Clamp(k,xvel_lo.z,xvel_hi.z)) <= zero) )
177  {
178  if ((dest_comp == RhoTheta_comp) && th_bc_ptr && th_file[3] != 0) {
179  // a theta profile reaches ext_dir only as rho*theta (th_file 1) or on an
180  // ext_dir_upwind face, which has no primitive variant (th_file 2):
181  // convert the latter with the ghost density here
182  const Real rho_th = (th_file[3] == 2)
183  ? (rho_is_extdir[3] ? rho_extdir_val[3] : dest_arr(dom_hi.x,j,k,Rho_comp))
184  : Real(1.0);
185  dest_arr(i,j,k,dest_comp) = rho_th * th_bc_ptr[k_profile];
186  } else {
187  dest_arr(i,j,k,dest_comp) = l_bc_extdir_vals_d[bc_comp][3];
188  }
189  } else if (h_bc_type == ERFBCType::ext_dir_prim) {
190  Real rho = rho_is_extdir[3] ? rho_extdir_val[3] : dest_arr(dom_hi.x,j,k,Rho_comp);
191  if ((dest_comp == RhoTheta_comp) && th_bc_ptr && th_file[3] == 2) {
192  dest_arr(i,j,k,dest_comp) = rho * th_bc_ptr[k_profile];
193  } else {
194  dest_arr(i,j,k,dest_comp) = rho * l_bc_extdir_vals_d[bc_comp][3];
195  }
196  }
197  }
198  );
199  }
200 
201  if (!is_periodic_in_y)
202  {
203  Real* th_bc_ptr = m_th_bc_data;
204  Box bx_ylo(bx); bx_ylo.setBig (1,dom_lo.y-1);
205  Box bx_yhi(bx); bx_yhi.setSmall(1,dom_hi.y+1);
206  //
207  // If we are setting Dirichlet values, set them in all ghost cells on an inflow face
208  // "bx" is already grown in the x- and y-directions so here we just grow it in z
209  //
210  bx_ylo.grow(2,ng[2]);
211  bx_yhi.grow(2,ng[2]);
212 
213  ParallelFor(
214  bx_ylo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n)
215  {
216  int dest_comp = icomp+n;
217  int k_profile = amrex::min(amrex::max(k, dom_lo.z), dom_hi.z) - dom_lo.z;
218  int bc_comp = (dest_comp >= RhoScalar_comp && dest_comp < RhoScalar_comp+NSCALARS) ?
219  BCVars::RhoScalar_bc_comp : dest_comp;
220  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
221  int l_bc_type = bc_ptr[n].lo(1);
222  if ( (l_bc_type == ERFBCType::ext_dir) ||
223  (l_bc_type == ERFBCType::ext_dir_upwind && yvel_arr(Clamp(i,yvel_lo.x,yvel_hi.x),dom_lo.y,Clamp(k,yvel_lo.z,yvel_hi.z)) >= zero) )
224  {
225  if ((dest_comp == RhoTheta_comp) && th_bc_ptr && th_file[1] != 0) {
226  // a theta profile reaches ext_dir only as rho*theta (th_file 1) or on an
227  // ext_dir_upwind face, which has no primitive variant (th_file 2):
228  // convert the latter with the ghost density here
229  const Real rho_th = (th_file[1] == 2)
230  ? (rho_is_extdir[1] ? rho_extdir_val[1] : dest_arr(i,dom_lo.y,k,Rho_comp))
231  : Real(1.0);
232  dest_arr(i,j,k,dest_comp) = rho_th * th_bc_ptr[k_profile];
233  } else {
234  dest_arr(i,j,k,dest_comp) = l_bc_extdir_vals_d[bc_comp][1];
235  }
236  } else if (l_bc_type == ERFBCType::ext_dir_prim) {
237  Real rho = rho_is_extdir[1] ? rho_extdir_val[1] : dest_arr(i,dom_lo.y,k,Rho_comp);
238  if ((dest_comp == RhoTheta_comp) && th_bc_ptr && th_file[1] == 2) {
239  dest_arr(i,j,k,dest_comp) = rho * th_bc_ptr[k_profile];
240  } else {
241  dest_arr(i,j,k,dest_comp) = rho * l_bc_extdir_vals_d[bc_comp][1];
242  }
243  }
244  },
245  bx_yhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n)
246  {
247  int dest_comp = icomp+n;
248  int k_profile = amrex::min(amrex::max(k, dom_lo.z), dom_hi.z) - dom_lo.z;
249  int bc_comp = (dest_comp >= RhoScalar_comp && dest_comp < RhoScalar_comp+NSCALARS) ?
250  BCVars::RhoScalar_bc_comp : dest_comp;
251  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
252  int h_bc_type = bc_ptr[n].hi(1);
253  if ( (h_bc_type == ERFBCType::ext_dir) ||
254  (h_bc_type == ERFBCType::ext_dir_upwind && yvel_arr(Clamp(i,yvel_lo.x,yvel_hi.x),dom_hi.y+1,Clamp(k,yvel_lo.z,yvel_hi.z)) <= zero) )
255  {
256  if ((dest_comp == RhoTheta_comp) && th_bc_ptr && th_file[4] != 0) {
257  // a theta profile reaches ext_dir only as rho*theta (th_file 1) or on an
258  // ext_dir_upwind face, which has no primitive variant (th_file 2):
259  // convert the latter with the ghost density here
260  const Real rho_th = (th_file[4] == 2)
261  ? (rho_is_extdir[4] ? rho_extdir_val[4] : dest_arr(i,dom_hi.y,k,Rho_comp))
262  : Real(1.0);
263  dest_arr(i,j,k,dest_comp) = rho_th * th_bc_ptr[k_profile];
264  } else {
265  dest_arr(i,j,k,dest_comp) = l_bc_extdir_vals_d[bc_comp][4];
266  }
267  } else if (h_bc_type == ERFBCType::ext_dir_prim) {
268  Real rho = rho_is_extdir[4] ? rho_extdir_val[4] : dest_arr(i,dom_hi.y,k,Rho_comp);
269  if ((dest_comp == RhoTheta_comp) && th_bc_ptr && th_file[4] == 2) {
270  dest_arr(i,j,k,dest_comp) = rho * th_bc_ptr[k_profile];
271  } else {
272  dest_arr(i,j,k,dest_comp) = rho * l_bc_extdir_vals_d[bc_comp][4];
273  }
274  }
275  }
276  );
277  }
278 
279  // Next do ghost cells in x-direction but not reaching out in y
280  // The corners we miss here will be covered in the y-loop below or by periodicity
281  if (!is_periodic_in_x)
282  {
283  // Populate ghost cells on lo-x and hi-x domain boundaries
284  Box bx_xlo(bx); bx_xlo.setBig (0,dom_lo.x-1);
285  Box bx_xhi(bx); bx_xhi.setSmall(0,dom_hi.x+1);
286  if (bx_xlo.smallEnd(2) != domain.smallEnd(2)) bx_xlo.growLo(2,ng[2]);
287  if (bx_xlo.bigEnd(2) != domain.bigEnd(2)) bx_xlo.growHi(2,ng[2]);
288  if (bx_xhi.smallEnd(2) != domain.smallEnd(2)) bx_xhi.growLo(2,ng[2]);
289  if (bx_xhi.bigEnd(2) != domain.bigEnd(2)) bx_xhi.growHi(2,ng[2]);
290  ParallelFor(
291  bx_xlo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n)
292  {
293  int dest_comp = icomp+n;
294  int l_bc_type = bc_ptr[n].lo(0);
295  int iflip = dom_lo.x - 1 - i;
296  if (l_bc_type == ERFBCType::foextrap) {
297  dest_arr(i,j,k,dest_comp) = dest_arr(dom_lo.x,j,k,dest_comp);
298  } else if (l_bc_type == ERFBCType::open) {
299  dest_arr(i,j,k,dest_comp) = dest_arr(dom_lo.x,j,k,dest_comp);
300  } else if (l_bc_type == ERFBCType::reflect_even) {
301  dest_arr(i,j,k,dest_comp) = dest_arr(iflip,j,k,dest_comp);
302  } else if (l_bc_type == ERFBCType::reflect_odd) {
303  dest_arr(i,j,k,dest_comp) = -dest_arr(iflip,j,k,dest_comp);
304  } else if (l_bc_type == ERFBCType::hoextrap) {
305  Real delta_i = static_cast<Real>(dom_lo.x - i);
306  dest_arr(i,j,k,dest_comp) = (one + delta_i)*dest_arr(dom_lo.x,j,k,dest_comp) - delta_i*dest_arr(dom_lo.x+1,j,k,dest_comp);
307  }
308  },
309  bx_xhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n)
310  {
311  int dest_comp = icomp+n;
312  int h_bc_type = bc_ptr[n].hi(0);
313  int iflip = 2*dom_hi.x + 1 - i;
314  if (h_bc_type == ERFBCType::foextrap) {
315  dest_arr(i,j,k,dest_comp) = dest_arr(dom_hi.x,j,k,dest_comp);
316  } else if (h_bc_type == ERFBCType::open) {
317  dest_arr(i,j,k,dest_comp) = dest_arr(dom_hi.x,j,k,dest_comp);
318  } else if (h_bc_type == ERFBCType::reflect_even) {
319  dest_arr(i,j,k,dest_comp) = dest_arr(iflip,j,k,dest_comp);
320  } else if (h_bc_type == ERFBCType::reflect_odd) {
321  dest_arr(i,j,k,dest_comp) = -dest_arr(iflip,j,k,dest_comp);
322  } else if (h_bc_type == ERFBCType::hoextrap) {
323  Real delta_i = static_cast<Real>(i - dom_hi.x);
324  dest_arr(i,j,k,dest_comp) = (one + delta_i)*dest_arr(dom_hi.x,j,k,dest_comp) - delta_i*dest_arr(dom_hi.x-1,j,k,dest_comp);
325  }
326  }
327  );
328  }
329 
330  if (!is_periodic_in_y)
331  {
332  // Populate ghost cells on lo-y and hi-y domain boundaries
333  Box bx_ylo(bx); bx_ylo.setBig (1,dom_lo.y-1);
334  Box bx_yhi(bx); bx_yhi.setSmall(1,dom_hi.y+1);
335  if (bx_ylo.smallEnd(2) != domain.smallEnd(2)) bx_ylo.growLo(2,ng[2]);
336  if (bx_ylo.bigEnd(2) != domain.bigEnd(2)) bx_ylo.growHi(2,ng[2]);
337  if (bx_yhi.smallEnd(2) != domain.smallEnd(2)) bx_yhi.growLo(2,ng[2]);
338  if (bx_yhi.bigEnd(2) != domain.bigEnd(2)) bx_yhi.growHi(2,ng[2]);
339  ParallelFor(
340  bx_ylo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n)
341  {
342  int dest_comp = icomp+n;
343  int l_bc_type = bc_ptr[n].lo(1);
344  int jflip = dom_lo.y - 1 - j;
345  if (l_bc_type == ERFBCType::foextrap) {
346  dest_arr(i,j,k,dest_comp) = dest_arr(i,dom_lo.y,k,dest_comp);
347  } else if (l_bc_type == ERFBCType::open) {
348  dest_arr(i,j,k,dest_comp) = dest_arr(i,dom_lo.y,k,dest_comp);
349  } else if (l_bc_type == ERFBCType::reflect_even) {
350  dest_arr(i,j,k,dest_comp) = dest_arr(i,jflip,k,dest_comp);
351  } else if (l_bc_type == ERFBCType::reflect_odd) {
352  dest_arr(i,j,k,dest_comp) = -dest_arr(i,jflip,k,dest_comp);
353  } else if (l_bc_type == ERFBCType::hoextrap) {
354  Real delta_j = static_cast<Real>(dom_lo.y - j);
355  dest_arr(i,j,k,dest_comp) = (one + delta_j)*dest_arr(i,dom_lo.y,k,dest_comp) - delta_j*dest_arr(i,dom_lo.y+1,k,dest_comp);
356  }
357 
358  },
359  bx_yhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n)
360  {
361  int dest_comp = icomp+n;
362  int h_bc_type = bc_ptr[n].hi(1);
363  int jflip = 2*dom_hi.y + 1 - j;
364  if (h_bc_type == ERFBCType::foextrap) {
365  dest_arr(i,j,k,dest_comp) = dest_arr(i,dom_hi.y,k,dest_comp);
366  } else if (h_bc_type == ERFBCType::open) {
367  dest_arr(i,j,k,dest_comp) = dest_arr(i,dom_hi.y,k,dest_comp);
368  } else if (h_bc_type == ERFBCType::reflect_even) {
369  dest_arr(i,j,k,dest_comp) = dest_arr(i,jflip,k,dest_comp);
370  } else if (h_bc_type == ERFBCType::reflect_odd) {
371  dest_arr(i,j,k,dest_comp) = -dest_arr(i,jflip,k,dest_comp);
372  } else if (h_bc_type == ERFBCType::hoextrap) {
373  Real delta_j = static_cast<Real>(j - dom_hi.y);
374  dest_arr(i,j,k,dest_comp) = (one + delta_j)*dest_arr(i,dom_hi.y,k,dest_comp) - delta_j*dest_arr(i,dom_hi.y-1,k,dest_comp);
375  }
376  }
377  );
378  }
379  Gpu::streamSynchronize();
380 }
OpenSide::hi
@ hi
OpenSide::lo
@ lo
RhoScalar_comp
#define RhoScalar_comp
Definition: ERF_IndexDefines.H:46
Rho_comp
#define Rho_comp
Definition: ERF_IndexDefines.H:42
RhoTheta_comp
#define RhoTheta_comp
Definition: ERF_IndexDefines.H:43
NBCVAR_max
#define NBCVAR_max
Definition: ERF_IndexDefines.H:32
NSCALARS
#define NSCALARS
Definition: ERF_IndexDefines.H:16
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);})
one
constexpr amrex::Real one
Definition: ERF_NumericalConstants.H:37
zero
constexpr amrex::Real zero
Definition: ERF_NumericalConstants.H:36
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:19
ERFPhysBCFunct_cons::impose_lateral_cons_bcs
void impose_lateral_cons_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 icomp, int ncomp, amrex::IntVect ng, const double time)
Definition: ERF_BoundaryConditionsCons.cpp:22
BCVars::RhoScalar_bc_comp
@ RhoScalar_bc_comp
Definition: ERF_IndexDefines.H:96
BCVars::Rho_bc_comp
@ Rho_bc_comp
Definition: ERF_IndexDefines.H:93
ERFBCType::open
@ open
Definition: ERF_IndexDefines.H:307
ERFBCType::reflect_odd
@ reflect_odd
Definition: ERF_IndexDefines.H:296
ERFBCType::hoextrap
@ hoextrap
Definition: ERF_IndexDefines.H:301
ERFBCType::foextrap
@ foextrap
Definition: ERF_IndexDefines.H:299
ERFBCType::ext_dir
@ ext_dir
Definition: ERF_IndexDefines.H:300
ERFBCType::ext_dir_prim
@ ext_dir_prim
Definition: ERF_IndexDefines.H:303
ERFBCType::ext_dir_upwind
@ ext_dir_upwind
Definition: ERF_IndexDefines.H:308
ERFBCType::int_dir
@ int_dir
Definition: ERF_IndexDefines.H:297
ERFBCType::reflect_even
@ reflect_even
Definition: ERF_IndexDefines.H:298
MicVar_Kess::rho
@ rho
Definition: ERF_Kessler.H:25
MicVar_Morr::ng
@ ng
Definition: ERF_Morrison.H:50
MicVar_Morr::nc
@ nc
Definition: ERF_Morrison.H:46
Here is the call graph for this function:

◆ impose_vertical_cons_bcs()

void ERFPhysBCFunct_cons::impose_vertical_cons_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  icomp,
int  ncomp,
const double  time,
bool  do_terrain_adjustment = true 
)

Impose vertical boundary conditions on conserved variables.

Parameters
[in,out]dest_arrcell-centered data to be filled
[in]bxbox holding data to be filled
[in]domainsimulation domain
[in]z_ndnodal physical height
[in]dxInvinverse cell size
[in]icompstarting component
[in]ncompnumber of components
[in]timetime at which data should be filled
[in]do_terrain_adjustmentwhether to apply terrain-aware Neumann adjustments

Impose vertical boundary conditions on conserved scalars (at cell centers)

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]icompthe index of the first component to be filled
[in]ncompthe number of components – this can be any value from 1 to NVAR as long as icomp+ncomp <= NVAR-one
[in]timetime at which the data should be filled
[in]do_terrain_adjustmentwhether to apply terrain-aware Neumann adjustments
403 {
404  BL_PROFILE_VAR("impose_vertical_cons_bcs()",impose_vertical_cons_bcs);
405  const auto& dom_lo = lbound(domain);
406  const auto& dom_hi = ubound(domain);
407 
408  Box per_grown_domain(domain);
409  int growx = (m_geom.isPeriodic(0)) ? 1 : 0;
410  int growy = (m_geom.isPeriodic(1)) ? 1 : 0;
411  per_grown_domain.grow(IntVect(growx,growy,0));
412  const auto& perdom_lo = lbound(per_grown_domain);
413  const auto& perdom_hi = ubound(per_grown_domain);
414 
415  GeometryData const& geomdata = m_geom.data();
416 
417  // xlo: ori = 0
418  // ylo: ori = 1
419  // zlo: ori = 2
420  // xhi: ori = 3
421  // yhi: ori = 4
422  // zhi: ori = 5
423 
424  // Based on BCRec for the domain, we need to make BCRec for this Box
425  // 0 is used as starting index for bcrs
426  Vector<BCRec> bcrs(ncomp);
427  GpuArray<GpuArray<Real, AMREX_SPACEDIM*2>,NBCVAR_max> l_bc_extdir_vals_d;
428  GpuArray<GpuArray<Real, AMREX_SPACEDIM*2>,NBCVAR_max> l_bc_neumann_vals_d;
429 
430  const int* bxlo = bx.loVect();
431  const int* bxhi = bx.hiVect();
432  const int* dlo = domain.loVect();
433  const int* dhi = domain.hiVect();
434 
435  for (int nc = 0; nc < ncomp; nc++)
436  {
437  int bc_comp = (icomp+nc >= RhoScalar_comp && icomp+nc < RhoScalar_comp+NSCALARS) ?
438  BCVars::RhoScalar_bc_comp : icomp+nc;
439  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
440  for (int dir = 0; dir < AMREX_SPACEDIM; dir++)
441  {
442  bcrs[nc].setLo(dir, ( bxlo[dir]<=dlo[dir]
443  ? m_domain_bcs_type[bc_comp].lo(dir) : BCType::int_dir ));
444  bcrs[nc].setHi(dir, ( bxhi[dir]>=dhi[dir]
445  ? m_domain_bcs_type[bc_comp].hi(dir) : BCType::int_dir ));
446  }
447 
448  for (int ori = 0; ori < 2*AMREX_SPACEDIM; ori++) {
449  l_bc_extdir_vals_d[bc_comp][ori] = m_bc_extdir_vals[bc_comp][ori];
450  l_bc_neumann_vals_d[bc_comp][ori] = m_bc_neumann_vals[bc_comp][ori];
451  }
452  }
453 
454  Gpu::DeviceVector<BCRec> bcrs_d(icomp+ncomp);
455  Gpu::copyAsync(Gpu::hostToDevice, bcrs.begin(), bcrs.end(), bcrs_d.begin());
456  const BCRec* bc_ptr = bcrs_d.data();
457 
458  {
459  Box bx_zlo(bx); bx_zlo.setBig (2,dom_lo.z-1);
460  Box bx_zhi(bx); bx_zhi.setSmall(2,dom_hi.z+1);
461  ParallelFor(
462  bx_zlo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n)
463  {
464  int dest_comp = icomp+n;
465  int bc_comp = (dest_comp >= RhoScalar_comp && dest_comp < RhoScalar_comp+NSCALARS) ?
466  BCVars::RhoScalar_bc_comp : dest_comp;
467  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
468  int l_bc_type = bc_ptr[n].lo(2);
469  if (l_bc_type == ERFBCType::ext_dir) {
470  dest_arr(i,j,k,dest_comp) = l_bc_extdir_vals_d[bc_comp][2];
471  } else if (l_bc_type == ERFBCType::ext_dir_prim) {
472  Real rho = dest_arr(i,j,dom_lo.z,Rho_comp);
473  dest_arr(i,j,k,dest_comp) = rho * l_bc_extdir_vals_d[bc_comp][2];
474  }
475  },
476  bx_zhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n)
477  {
478  int dest_comp = icomp+n;
479  int bc_comp = (dest_comp >= RhoScalar_comp && dest_comp < RhoScalar_comp+NSCALARS) ?
480  BCVars::RhoScalar_bc_comp : dest_comp;
481  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
482  int h_bc_type = bc_ptr[n].hi(2);
483  if (h_bc_type == ERFBCType::ext_dir) {
484  dest_arr(i,j,k,dest_comp) = l_bc_extdir_vals_d[bc_comp][5];
485  } else if (h_bc_type == ERFBCType::ext_dir_prim) {
486  Real rho = dest_arr(i,j,dom_hi.z,Rho_comp);
487  dest_arr(i,j,k,dest_comp) = rho * l_bc_extdir_vals_d[bc_comp][5];
488  }
489 
490  }
491  );
492  }
493 
494  {
495  Box bx_zlo(bx); bx_zlo.setBig (2,dom_lo.z-1);
496  Box bx_zhi(bx); bx_zhi.setSmall(2,dom_hi.z+1);
497  // Populate ghost cells on lo-z and hi-z domain boundaries
498  ParallelFor(
499  bx_zlo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n)
500  {
501  int dest_comp = icomp+n;
502  int bc_comp = (dest_comp >= RhoScalar_comp && dest_comp < RhoScalar_comp+NSCALARS) ?
503  BCVars::RhoScalar_bc_comp : dest_comp;
504  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
505  int l_bc_type = bc_ptr[n].lo(2);
506  int kflip = dom_lo.z - 1 - k;
507  if (l_bc_type == ERFBCType::foextrap) {
508  dest_arr(i,j,k,dest_comp) = dest_arr(i,j,dom_lo.z,dest_comp);
509  } else if (l_bc_type == ERFBCType::open) {
510  dest_arr(i,j,k,dest_comp) = dest_arr(i,j,dom_lo.z,dest_comp);
511  } else if (l_bc_type == ERFBCType::reflect_even) {
512  dest_arr(i,j,k,dest_comp) = dest_arr(i,j,kflip,dest_comp);
513  } else if (l_bc_type == ERFBCType::reflect_odd) {
514  dest_arr(i,j,k,dest_comp) = -dest_arr(i,j,kflip,dest_comp);
515  } else if (l_bc_type == ERFBCType::neumann) {
516  Real delta_z = Compute_Z_AtCellCenter(i,j,dom_lo.z,z_phys_nd)
517  - Compute_Z_AtCellCenter(i,j,k ,z_phys_nd);
518  if( (icomp+n) == Rho_comp ) {
519  dest_arr(i,j,k,dest_comp) = dest_arr(i,j,dom_lo.z,dest_comp) -
520  delta_z*l_bc_neumann_vals_d[bc_comp][2];
521  } else {
522  dest_arr(i,j,k,dest_comp) = dest_arr(i,j,dom_lo.z,dest_comp) -
523  delta_z*l_bc_neumann_vals_d[bc_comp][2]*dest_arr(i,j,dom_lo.z,Rho_comp);
524  }
525  } else if (l_bc_type == ERFBCType::hoextrap) {
526  Real delta_k = static_cast<Real>(dom_lo.z - k);
527  dest_arr(i,j,k,dest_comp) = (one + delta_k) * dest_arr(i,j,dom_lo.z ,dest_comp) -
528  delta_k * dest_arr(i,j,dom_lo.z+1,dest_comp);
529  }
530  },
531  bx_zhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n)
532  {
533  int dest_comp = icomp+n;
534  int bc_comp = (dest_comp >= RhoScalar_comp && dest_comp < RhoScalar_comp+NSCALARS) ?
535  BCVars::RhoScalar_bc_comp : dest_comp;
536  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
537  int h_bc_type = bc_ptr[n].hi(2);
538  int kflip = 2*dom_hi.z + 1 - k;
539  if (h_bc_type == ERFBCType::foextrap) {
540  dest_arr(i,j,k,dest_comp) = dest_arr(i,j,dom_hi.z,dest_comp);
541  } else if (h_bc_type == ERFBCType::open) {
542  dest_arr(i,j,k,dest_comp) = dest_arr(i,j,dom_hi.z,dest_comp);
543  } else if (h_bc_type == ERFBCType::reflect_even) {
544  dest_arr(i,j,k,dest_comp) = dest_arr(i,j,kflip,dest_comp);
545  } else if (h_bc_type == ERFBCType::reflect_odd) {
546  dest_arr(i,j,k,dest_comp) = -dest_arr(i,j,kflip,dest_comp);
547  } else if (h_bc_type == ERFBCType::neumann) {
548  Real delta_z = Compute_Z_AtCellCenter(i,j,k ,z_phys_nd)
549  - Compute_Z_AtCellCenter(i,j,dom_hi.z,z_phys_nd);
550  if( (icomp+n) == Rho_comp ) {
551  dest_arr(i,j,k,dest_comp) = dest_arr(i,j,dom_hi.z,dest_comp) +
552  delta_z*l_bc_neumann_vals_d[bc_comp][5];
553  } else {
554  dest_arr(i,j,k,dest_comp) = dest_arr(i,j,dom_hi.z,dest_comp) +
555  delta_z*l_bc_neumann_vals_d[bc_comp][5]*dest_arr(i,j,dom_hi.z,Rho_comp);
556  }
557  } else if (h_bc_type == ERFBCType::hoextrap){
558  Real delta_k = static_cast<Real>(k - dom_hi.z);
559  dest_arr(i,j,k,dest_comp) = (one + delta_k)*dest_arr(i,j,dom_hi.z,dest_comp) - delta_k*dest_arr(i,j,dom_hi.z-1,dest_comp);
560  }
561  }
562  );
563  }
564 
565  if (do_terrain_adjustment && m_z_phys_nd) {
566  const auto& bx_lo = lbound(bx);
567  const auto& bx_hi = ubound(bx);
568  const BCRec* bc_ptr_h = bcrs.data();
569  // Neumann conditions (d<var>/dn = 0) must be aware of the surface normal with terrain.
570  // An additional source term arises from d<var>/dx & d<var>/dy & met_h_xi/eta/zeta.
571  //=====================================================================================
572  // Only modify scalars, U, or V
573  // Loop over each component
574  for (int n = 0; n < ncomp; n++) {
575  // Hit for Neumann condition at kmin
576  int dest_comp = icomp+n;
577  int l_bc_type = bc_ptr_h[n].lo(2);
578  if(l_bc_type == ERFBCType::foextrap)
579  {
580  // Loop over ghost cells in bottom XY-plane (valid box)
581  Box xybx = bx;
582 
583  int k0 = dom_lo.z;
584  if (xybx.smallEnd(2) < 0)
585  {
586  xybx.setBig(2,dom_lo.z-1);
587  xybx.setSmall(2,bx.smallEnd(2));
588 
589  // Get the dz cell size
590  Real dz = geomdata.CellSize(2);
591 
592  // Fill all the Neumann srcs with terrain
593  ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
594  {
595  // Clip indices for ghost-cells
596  int ii = amrex::min(amrex::max(i,perdom_lo.x),perdom_hi.x);
597  int jj = amrex::min(amrex::max(j,perdom_lo.y),perdom_hi.y);
598 
599  // Get metrics
600  Real met_h_xi = Compute_h_xi_AtCellCenter (ii,jj,k0,dxInv,z_phys_nd);
601  Real met_h_eta = Compute_h_eta_AtCellCenter (ii,jj,k0,dxInv,z_phys_nd);
602  Real met_h_zeta = Compute_h_zeta_AtCellCenter(ii,jj,k0,dxInv,z_phys_nd);
603 
604  // GradX at IJK location inside domain -- this relies on the assumption that we have
605  // used foextrap for cell-centered quantities outside the domain to define the gradient as zero
606  Real GradVarx, GradVary;
607  if (i < dom_lo.x-1 || i > dom_hi.x+1 || (i+1 > bx_hi.x && i-1 < bx_lo.x) ) {
608  GradVarx = zero;
609  } else if (i+1 > bx_hi.x) {
610  GradVarx = dxInv[0] * (dest_arr(i ,j,k0,dest_comp) - dest_arr(i-1,j,k0,dest_comp));
611  } else if (i-1 < bx_lo.x) {
612  GradVarx = dxInv[0] * (dest_arr(i+1,j,k0,dest_comp) - dest_arr(i ,j,k0,dest_comp));
613  } else {
614  GradVarx = myhalf * dxInv[0] * (dest_arr(i+1,j,k0,dest_comp) - dest_arr(i-1,j,k0,dest_comp));
615  }
616 
617  // GradY at IJK location inside domain -- this relies on the assumption that we have
618  // used foextrap for cell-centered quantities outside the domain to define the gradient as zero
619  if (j < dom_lo.y-1 || j > dom_hi.y+1 || (j+1 > bx_hi.y && j-1 < bx_lo.y) ) {
620  GradVary = zero;
621  } else if (j+1 > bx_hi.y) {
622  GradVary = dxInv[1] * (dest_arr(i,j ,k0,dest_comp) - dest_arr(i,j-1,k0,dest_comp));
623  } else if (j-1 < bx_lo.y) {
624  GradVary = dxInv[1] * (dest_arr(i,j+1,k0,dest_comp) - dest_arr(i,j ,k0,dest_comp));
625  } else {
626  GradVary = myhalf * dxInv[1] * (dest_arr(i,j+1,k0,dest_comp) - dest_arr(i,j-1,k0,dest_comp));
627  }
628 
629  // Prefactor
630  Real met_fac = met_h_zeta / ( met_h_xi*met_h_xi + met_h_eta*met_h_eta + one );
631 
632  // Accumulate in bottom ghost cell (EXTRAP already populated)
633  dest_arr(i,j,k,dest_comp) -= dz * met_fac * ( met_h_xi * GradVarx + met_h_eta * GradVary );
634  });
635  } // box includes k0
636  } // foextrap
637  } // ncomp
638  } // m_z_phys_nd
639  Gpu::streamSynchronize();
640 }
dxInv
amrex::GpuArray< Real, AMREX_SPACEDIM > dxInv
Definition: ERF_InitCustomPertVels_ParticleTests.H:17
myhalf
constexpr amrex::Real myhalf
Definition: ERF_NumericalConstants.H:41
Compute_Z_AtCellCenter
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_Z_AtCellCenter(const int &i, const int &j, const int &k, const amrex::Array4< const amrex::Real > &z_nd)
Definition: ERF_TerrainMetrics.H:712
Compute_h_xi_AtCellCenter
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_xi_AtCellCenter(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:222
Compute_h_eta_AtCellCenter
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_eta_AtCellCenter(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:247
Compute_h_zeta_AtCellCenter
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_zeta_AtCellCenter(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:197
ERFPhysBCFunct_cons::impose_vertical_cons_bcs
void impose_vertical_cons_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 icomp, int ncomp, const double time, bool do_terrain_adjustment=true)
Definition: ERF_BoundaryConditionsCons.cpp:397
ERFBCType::neumann
@ neumann
Definition: ERF_IndexDefines.H:305
WDM6SedCellScratch::dz
@ dz
Definition: ERF_AdvanceWDM6.cpp:272
Here is the call graph for this function:

◆ operator()()

void ERFPhysBCFunct_cons::operator() ( amrex::MultiFab &  mf,
amrex::MultiFab &  xvel,
amrex::MultiFab &  yvel,
int  icomp,
int  ncomp,
amrex::IntVect const &  nghost,
const double  time,
int  bccomp_cons,
bool  do_fb = true,
bool  do_terrain_adjustment = true 
)

Impose physical boundary conditions at domain boundaries

Parameters
[in,out]mfconserved-variable MultiFab to be filled
[in]xvelx-velocity MultiFab used by upwind boundary conditions
[in]yvely-velocity MultiFab used by upwind boundary conditions
[in]icompstarting component for conserved variables
[in]ncompnumber of components for conserved variables
[in]nghostnumber of ghost cells to be filled for conserved variables
[in]timetime at which the data should be filled
[in]bccomp_consindex into m_domain_bcs_type for conserved variables
[in]do_fbwhether to fill periodic/interior ghost cells first
[in]do_terrain_adjustmentwhether to apply terrain-aware Neumann adjustments

Impose physical boundary conditions at domain boundaries

Parameters
[in,out]mfMultiFab of cell-centered quantities to be filled
[in]xvelx-velocity MultiFab used by upwind boundary conditions
[in]yvely-velocity MultiFab used by upwind boundary conditions
[in]icompstarting component
[in]ncompnumber of components
[in]nghostnumber of ghost cells to be filled
[in]timetime at which data should be filled
[in]bccompindex into m_domain_bcs_type, unused here
[in]do_fbwhether to fill periodic/interior ghost cells first
[in]do_terrain_adjustmentwhether to apply terrain-aware Neumann adjustments
91 {
92  BL_PROFILE("ERFPhysBCFunct_cons::()");
93 
94  //
95  // We fill all of the interior and periodic ghost cells first, so we can fill
96  // those directly inside the lateral and vertical calls.
97  // If triply periodic this is all we do
98  //
99  if (do_fb) {
100  mf.FillBoundary(icomp,ncomp,m_geom.periodicity());
101  }
102 
103  if (m_geom.isAllPeriodic()) return;
104 
105  const auto& domain = m_geom.Domain();
106  const auto dxInv = m_geom.InvCellSizeArray();
107 
108  // Create a grown domain box containing valid + periodic cells
109  Box gdomain = domain;
110  for (int i = 0; i < AMREX_SPACEDIM; ++i) {
111  if (m_geom.isPeriodic(i)) {
112  gdomain.grow(i, nghost[i]);
113  }
114  }
115 
116 #ifdef AMREX_USE_OMP
117 #pragma omp parallel if (Gpu::notInLaunchRegion())
118 #endif
119  {
120  for (MFIter mfi(mf,false); mfi.isValid(); ++mfi)
121  {
122  //
123  // This is the box we pass to the different routines
124  // NOTE -- this is the full grid box NOT the tile box
125  //
126  Box bx = mfi.validbox();
127 
128  //
129  // These are the boxes we use to test on relative to the domain
130  //
131  Box cbx1 = bx; cbx1.grow(IntVect(nghost[0],nghost[1],0));
132  Box cbx2 = bx; cbx2.grow(nghost);
133 
134  Array4<const Real> z_nd_arr;
135 
136  if (m_z_phys_nd) {
137  z_nd_arr = m_z_phys_nd->const_array(mfi);
138  }
139 
140  if (!gdomain.contains(cbx2)) {
141  Array4< Real> const& cons_arr = mf.array(mfi);
142  Array4<const Real> const& velx_arr = xvel.const_array(mfi);
143  Array4<const Real> const& vely_arr = yvel.const_array(mfi);
144 
145  if (!m_use_real_bcs) {
146  // We send a box with ghost cells in the lateral directions only
147  impose_lateral_cons_bcs(cons_arr,velx_arr,vely_arr,cbx1,domain,icomp,ncomp,nghost,time);
148  }
149 
150  // We send the full FAB box with ghost cells
151  impose_vertical_cons_bcs(cons_arr,cbx2,domain,z_nd_arr,dxInv,icomp,ncomp,time,do_terrain_adjustment);
152  }
153 
154  } // MFIter
155  } // OpenMP
156 
157  // The terrain correction below the ground is computed box by box; make the
158  // copies of each ghost cell agree. Only a terrain-fitted mesh has a
159  // non-zero correction -- z_phys_nd is allocated (and zeroed) for every
160  // terrain type, so testing it alone would pay for the sync on flat runs.
161  if (do_terrain_adjustment && m_z_phys_nd && terrain_is_fitted(m_terrain_type)) {
162  m_below_ground_sync(mf, icomp, ncomp, nghost, m_geom);
163  }
164 } // operator()
ERFPhysBCFunct_cons::m_below_ground_sync
BelowGroundGhostSync m_below_ground_sync
Definition: ERF_PhysBCFunct.H:171
Vars::xvel
@ xvel
Definition: ERF_IndexDefines.H:218
Vars::yvel
@ yvel
Definition: ERF_IndexDefines.H:219

Member Data Documentation

◆ m_bc_extdir_vals

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

◆ m_bc_neumann_vals

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

◆ m_below_ground_sync

BelowGroundGhostSync ERFPhysBCFunct_cons::m_below_ground_sync
private

◆ m_domain_bcs_type

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

◆ m_domain_bcs_type_d

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

◆ m_geom

amrex::Geometry ERFPhysBCFunct_cons::m_geom
private

◆ m_lev

int ERFPhysBCFunct_cons::m_lev
private

◆ m_terrain_type

TerrainType ERFPhysBCFunct_cons::m_terrain_type
private

◆ m_th_bc_data

amrex::Real* ERFPhysBCFunct_cons::m_th_bc_data
private

◆ m_th_file_face

amrex::GpuArray<int,AMREX_SPACEDIM*2> ERFPhysBCFunct_cons::m_th_file_face
private

Referenced by ERFPhysBCFunct_cons().

◆ m_use_real_bcs

bool ERFPhysBCFunct_cons::m_use_real_bcs
private

◆ m_z_phys_nd

amrex::MultiFab* ERFPhysBCFunct_cons::m_z_phys_nd
private
The documentation for this class was generated from the following files: