ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
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, std::unique_ptr< amrex::MultiFab > &z_phys_nd, const bool use_real_bcs, amrex::Real *th_bc_data)
 
 ~ERFPhysBCFunct_cons ()
 
void operator() (amrex::MultiFab &mf, amrex::MultiFab &xvel, amrex::MultiFab &yvel, int icomp, int ncomp, amrex::IntVect const &nghost, const amrex::Real 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)
 
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, 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
 
amrex::MultiFab * m_z_phys_nd
 
bool m_use_real_bcs
 
amrex::Real * m_th_bc_data
 

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,
std::unique_ptr< amrex::MultiFab > &  z_phys_nd,
const bool  use_real_bcs,
amrex::Real *  th_bc_data 
)
inline
34  : m_lev(lev), m_geom(geom),
35  m_domain_bcs_type(domain_bcs_type),
36  m_domain_bcs_type_d(domain_bcs_type_d),
37  m_bc_extdir_vals(bc_extdir_vals),
38  m_bc_neumann_vals(bc_neumann_vals),
39  m_z_phys_nd(z_phys_nd.get()),
40  m_use_real_bcs(use_real_bcs),
41  m_th_bc_data(th_bc_data)
42  {}
ERFPhysBCFunct_cons::m_use_real_bcs
bool m_use_real_bcs
Definition: ERF_PhysBCFunct.H:79
ERFPhysBCFunct_cons::m_geom
amrex::Geometry m_geom
Definition: ERF_PhysBCFunct.H:73
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:76
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:77
ERFPhysBCFunct_cons::m_domain_bcs_type_d
amrex::Gpu::DeviceVector< amrex::BCRec > m_domain_bcs_type_d
Definition: ERF_PhysBCFunct.H:75
ERFPhysBCFunct_cons::m_domain_bcs_type
amrex::Vector< amrex::BCRec > m_domain_bcs_type
Definition: ERF_PhysBCFunct.H:74
ERFPhysBCFunct_cons::m_th_bc_data
amrex::Real * m_th_bc_data
Definition: ERF_PhysBCFunct.H:80
ERFPhysBCFunct_cons::m_z_phys_nd
amrex::MultiFab * m_z_phys_nd
Definition: ERF_PhysBCFunct.H:78
ERFPhysBCFunct_cons::m_lev
int m_lev
Definition: ERF_PhysBCFunct.H:72

◆ ~ERFPhysBCFunct_cons()

ERFPhysBCFunct_cons::~ERFPhysBCFunct_cons ( )
inline
44 {}

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 
)
22 {
23  BL_PROFILE_VAR("impose_lateral_cons_bcs()",impose_lateral_cons_bcs);
24  const auto& dom_lo = lbound(domain);
25  const auto& dom_hi = ubound(domain);
26 
27  // xlo: ori = 0
28  // ylo: ori = 1
29  // zlo: ori = 2
30  // xhi: ori = 3
31  // yhi: ori = 4
32  // zhi: ori = 5
33 
34  // Based on BCRec for the domain, we need to make BCRec for this Box
35  // 0 is used as starting index for bcrs
36  Vector<BCRec> bcrs(ncomp);
37 
38  GpuArray<GpuArray<Real, AMREX_SPACEDIM*2>,NBCVAR_max> l_bc_extdir_vals_d;
39 
40  const int* bxlo = bx.loVect();
41  const int* bxhi = bx.hiVect();
42  const int* dlo = domain.loVect();
43  const int* dhi = domain.hiVect();
44 
45  for (int nc = 0; nc < ncomp; nc++)
46  {
47  int bc_comp = (icomp+nc >= RhoScalar_comp && icomp+nc < RhoScalar_comp+NSCALARS) ?
48  BCVars::RhoScalar_bc_comp : icomp+nc;
49  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
50  for (int dir = 0; dir < AMREX_SPACEDIM; dir++)
51  {
52  bcrs[nc].setLo(dir, ( bxlo[dir]<=dlo[dir]
53  ? m_domain_bcs_type[bc_comp].lo(dir) : BCType::int_dir ));
54  bcrs[nc].setHi(dir, ( bxhi[dir]>=dhi[dir]
55  ? m_domain_bcs_type[bc_comp].hi(dir) : BCType::int_dir ));
56  }
57 
58  for (int ori = 0; ori < 2*AMREX_SPACEDIM; ori++) {
59  l_bc_extdir_vals_d[bc_comp][ori] = m_bc_extdir_vals[bc_comp][ori];
60  }
61  }
62 
63  Gpu::DeviceVector<BCRec> bcrs_d(ncomp);
64  Gpu::copyAsync(Gpu::hostToDevice, bcrs.begin(), bcrs.end(), bcrs_d.begin());
65  const BCRec* bc_ptr = bcrs_d.data();
66 
67  GeometryData const& geomdata = m_geom.data();
68  bool is_periodic_in_x = geomdata.isPeriodic(0);
69  bool is_periodic_in_y = geomdata.isPeriodic(1);
70 
71  // First do all ext_dir bcs
72  if (!is_periodic_in_x)
73  {
74  Real* th_bc_ptr = m_th_bc_data;
75  Box bx_xlo(bx); bx_xlo.setBig (0,dom_lo.x-1);
76  Box bx_xhi(bx); bx_xhi.setSmall(0,dom_hi.x+1);
77  //
78  // If we are setting Dirichlet values, set them in all ghost cells on an inflow face
79  // "bx" is already grown in the x- and y-directions so here we just grow it in z
80  //
81  bx_xlo.grow(2,ng[2]);
82  bx_xhi.grow(2,ng[2]);
83  ParallelFor(
84  bx_xlo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n)
85  {
86  int dest_comp = icomp+n;
87  int bc_comp = (dest_comp >= RhoScalar_comp && dest_comp < RhoScalar_comp+NSCALARS) ?
88  BCVars::RhoScalar_bc_comp : dest_comp;
89  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
90  int l_bc_type = bc_ptr[n].lo(0);
91 
92  if ( (l_bc_type == ERFBCType::ext_dir) ||
93  (l_bc_type == ERFBCType::ext_dir_upwind && xvel_arr(dom_lo.x,j,k) >= 0.) )
94  {
95  if ((dest_comp == RhoTheta_comp) && th_bc_ptr) {
96  dest_arr(i,j,k,dest_comp) = th_bc_ptr[k];
97  } else {
98  dest_arr(i,j,k,dest_comp) = l_bc_extdir_vals_d[bc_comp][0];
99  }
100  } else if (l_bc_type == ERFBCType::ext_dir_prim) {
101  Real rho = dest_arr(dom_lo.x,j,k,Rho_comp);
102  if ((dest_comp == RhoTheta_comp) && th_bc_ptr) {
103  dest_arr(i,j,k,dest_comp) = rho * th_bc_ptr[k];
104  } else {
105  dest_arr(i,j,k,dest_comp) = rho * l_bc_extdir_vals_d[bc_comp][0];
106  }
107  }
108  },
109  bx_xhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n)
110  {
111  int dest_comp = icomp+n;
112  int bc_comp = (dest_comp >= RhoScalar_comp && dest_comp < RhoScalar_comp+NSCALARS) ?
113  BCVars::RhoScalar_bc_comp : dest_comp;
114  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
115  int h_bc_type = bc_ptr[n].hi(0);
116 
117  if ( (h_bc_type == ERFBCType::ext_dir) ||
118  (h_bc_type == ERFBCType::ext_dir_upwind && xvel_arr(dom_hi.x+1,j,k) <= 0.) )
119  {
120  if ((dest_comp == RhoTheta_comp) && th_bc_ptr) {
121  dest_arr(i,j,k,dest_comp) = th_bc_ptr[k];
122  } else {
123  dest_arr(i,j,k,dest_comp) = l_bc_extdir_vals_d[bc_comp][3];
124  }
125  } else if (h_bc_type == ERFBCType::ext_dir_prim) {
126  Real rho = dest_arr(dom_hi.x,j,k,Rho_comp);
127  if ((dest_comp == RhoTheta_comp) && th_bc_ptr) {
128  dest_arr(i,j,k,dest_comp) = rho * th_bc_ptr[k];
129  } else {
130  dest_arr(i,j,k,dest_comp) = rho * l_bc_extdir_vals_d[bc_comp][3];
131  }
132  }
133  }
134  );
135  }
136 
137  if (!is_periodic_in_y)
138  {
139  Real* th_bc_ptr = m_th_bc_data;
140  Box bx_ylo(bx); bx_ylo.setBig (1,dom_lo.y-1);
141  Box bx_yhi(bx); bx_yhi.setSmall(1,dom_hi.y+1);
142  //
143  // If we are setting Dirichlet values, set them in all ghost cells on an inflow face
144  // "bx" is already grown in the x- and y-directions so here we just grow it in z
145  //
146  bx_ylo.grow(2,ng[2]);
147  bx_yhi.grow(2,ng[2]);
148 
149  ParallelFor(
150  bx_ylo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n)
151  {
152  int dest_comp = icomp+n;
153  int bc_comp = (dest_comp >= RhoScalar_comp && dest_comp < RhoScalar_comp+NSCALARS) ?
154  BCVars::RhoScalar_bc_comp : dest_comp;
155  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
156  int l_bc_type = bc_ptr[n].lo(1);
157  if ( (l_bc_type == ERFBCType::ext_dir) ||
158  (l_bc_type == ERFBCType::ext_dir_upwind && yvel_arr(i,dom_lo.y,k) >= 0.) )
159  {
160  if ((dest_comp == RhoTheta_comp) && th_bc_ptr) {
161  dest_arr(i,j,k,dest_comp) = th_bc_ptr[k];
162  } else {
163  dest_arr(i,j,k,dest_comp) = l_bc_extdir_vals_d[bc_comp][1];
164  }
165  } else if (l_bc_type == ERFBCType::ext_dir_prim) {
166  Real rho = dest_arr(i,dom_lo.y,k,Rho_comp);
167  if ((dest_comp == RhoTheta_comp) && th_bc_ptr) {
168  dest_arr(i,j,k,dest_comp) = rho * th_bc_ptr[k];
169  } else {
170  dest_arr(i,j,k,dest_comp) = rho * l_bc_extdir_vals_d[bc_comp][1];
171  }
172  }
173  },
174  bx_yhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n)
175  {
176  int dest_comp = icomp+n;
177  int bc_comp = (dest_comp >= RhoScalar_comp && dest_comp < RhoScalar_comp+NSCALARS) ?
178  BCVars::RhoScalar_bc_comp : dest_comp;
179  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
180  int h_bc_type = bc_ptr[n].hi(1);
181  if ( (h_bc_type == ERFBCType::ext_dir) ||
182  (h_bc_type == ERFBCType::ext_dir_upwind && yvel_arr(i,dom_hi.y+1,k) <= 0.) )
183  {
184  if ((dest_comp == RhoTheta_comp) && th_bc_ptr) {
185  dest_arr(i,j,k,dest_comp) = th_bc_ptr[k];
186  } else {
187  dest_arr(i,j,k,dest_comp) = l_bc_extdir_vals_d[bc_comp][4];
188  }
189  } else if (h_bc_type == ERFBCType::ext_dir_prim) {
190  Real rho = dest_arr(i,dom_hi.y,k,Rho_comp);
191  if ((dest_comp == RhoTheta_comp) && th_bc_ptr) {
192  dest_arr(i,j,k,dest_comp) = rho * th_bc_ptr[k];
193  } else {
194  dest_arr(i,j,k,dest_comp) = rho * l_bc_extdir_vals_d[bc_comp][4];
195  }
196  }
197  }
198  );
199  }
200 
201  // Next do ghost cells in x-direction but not reaching out in y
202  // The corners we miss here will be covered in the y-loop below or by periodicity
203  if (!is_periodic_in_x)
204  {
205  // Populate ghost cells on lo-x and hi-x domain boundaries
206  Box bx_xlo(bx); bx_xlo.setBig (0,dom_lo.x-1);
207  Box bx_xhi(bx); bx_xhi.setSmall(0,dom_hi.x+1);
208  if (bx_xlo.smallEnd(2) != domain.smallEnd(2)) bx_xlo.growLo(2,ng[2]);
209  if (bx_xlo.bigEnd(2) != domain.bigEnd(2)) bx_xlo.growHi(2,ng[2]);
210  if (bx_xhi.smallEnd(2) != domain.smallEnd(2)) bx_xhi.growLo(2,ng[2]);
211  if (bx_xhi.bigEnd(2) != domain.bigEnd(2)) bx_xhi.growHi(2,ng[2]);
212  ParallelFor(
213  bx_xlo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n)
214  {
215  int dest_comp = icomp+n;
216  int l_bc_type = bc_ptr[n].lo(0);
217  int iflip = dom_lo.x - 1 - i;
218  if (l_bc_type == ERFBCType::foextrap) {
219  dest_arr(i,j,k,dest_comp) = dest_arr(dom_lo.x,j,k,dest_comp);
220  } else if (l_bc_type == ERFBCType::open) {
221  dest_arr(i,j,k,dest_comp) = dest_arr(dom_lo.x,j,k,dest_comp);
222  } else if (l_bc_type == ERFBCType::reflect_even) {
223  dest_arr(i,j,k,dest_comp) = dest_arr(iflip,j,k,dest_comp);
224  } else if (l_bc_type == ERFBCType::reflect_odd) {
225  dest_arr(i,j,k,dest_comp) = -dest_arr(iflip,j,k,dest_comp);
226  } else if (l_bc_type == ERFBCType::hoextrapcc) {
227  Real delta_i = (dom_lo.x - i);
228  dest_arr(i,j,k,dest_comp) = (1.0 + delta_i)*dest_arr(dom_lo.x,j,k,dest_comp) - delta_i*dest_arr(dom_lo.x+1,j,k,dest_comp) ;
229  }
230  },
231  bx_xhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n)
232  {
233  int dest_comp = icomp+n;
234  int h_bc_type = bc_ptr[n].hi(0);
235  int iflip = 2*dom_hi.x + 1 - i;
236  if (h_bc_type == ERFBCType::foextrap) {
237  dest_arr(i,j,k,dest_comp) = dest_arr(dom_hi.x,j,k,dest_comp);
238  } else if (h_bc_type == ERFBCType::open) {
239  dest_arr(i,j,k,dest_comp) = dest_arr(dom_hi.x,j,k,dest_comp);
240  } else if (h_bc_type == ERFBCType::reflect_even) {
241  dest_arr(i,j,k,dest_comp) = dest_arr(iflip,j,k,dest_comp);
242  } else if (h_bc_type == ERFBCType::reflect_odd) {
243  dest_arr(i,j,k,dest_comp) = -dest_arr(iflip,j,k,dest_comp);
244  } else if (h_bc_type == ERFBCType::hoextrapcc) {
245  Real delta_i = (i - dom_hi.x);
246  dest_arr(i,j,k,dest_comp) = (1.0 + delta_i)*dest_arr(dom_hi.x,j,k,dest_comp) - delta_i*dest_arr(dom_hi.x-1,j,k,dest_comp) ;
247  }
248  }
249  );
250  }
251 
252  if (!is_periodic_in_y)
253  {
254  // Populate ghost cells on lo-y and hi-y domain boundaries
255  Box bx_ylo(bx); bx_ylo.setBig (1,dom_lo.y-1);
256  Box bx_yhi(bx); bx_yhi.setSmall(1,dom_hi.y+1);
257  if (bx_ylo.smallEnd(2) != domain.smallEnd(2)) bx_ylo.growLo(2,ng[2]);
258  if (bx_ylo.bigEnd(2) != domain.bigEnd(2)) bx_ylo.growHi(2,ng[2]);
259  if (bx_yhi.smallEnd(2) != domain.smallEnd(2)) bx_yhi.growLo(2,ng[2]);
260  if (bx_yhi.bigEnd(2) != domain.bigEnd(2)) bx_yhi.growHi(2,ng[2]);
261  ParallelFor(
262  bx_ylo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n)
263  {
264  int dest_comp = icomp+n;
265  int l_bc_type = bc_ptr[n].lo(1);
266  int jflip = dom_lo.y - 1 - j;
267  if (l_bc_type == ERFBCType::foextrap) {
268  dest_arr(i,j,k,dest_comp) = dest_arr(i,dom_lo.y,k,dest_comp);
269  } else if (l_bc_type == ERFBCType::open) {
270  dest_arr(i,j,k,dest_comp) = dest_arr(i,dom_lo.y,k,dest_comp);
271  } else if (l_bc_type == ERFBCType::reflect_even) {
272  dest_arr(i,j,k,dest_comp) = dest_arr(i,jflip,k,dest_comp);
273  } else if (l_bc_type == ERFBCType::reflect_odd) {
274  dest_arr(i,j,k,dest_comp) = -dest_arr(i,jflip,k,dest_comp);
275  } else if (l_bc_type == ERFBCType::hoextrapcc) {
276  Real delta_j = (dom_lo.y - j);
277  dest_arr(i,j,k,dest_comp) = (1.0 + delta_j)*dest_arr(i,dom_lo.y,k,dest_comp) - delta_j*dest_arr(i,dom_lo.y+1,k,dest_comp) ;
278  }
279 
280  },
281  bx_yhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n)
282  {
283  int dest_comp = icomp+n;
284  int h_bc_type = bc_ptr[n].hi(1);
285  int jflip = 2*dom_hi.y + 1 - j;
286  if (h_bc_type == ERFBCType::foextrap) {
287  dest_arr(i,j,k,dest_comp) = dest_arr(i,dom_hi.y,k,dest_comp);
288  } else if (h_bc_type == ERFBCType::open) {
289  dest_arr(i,j,k,dest_comp) = dest_arr(i,dom_hi.y,k,dest_comp);
290  } else if (h_bc_type == ERFBCType::reflect_even) {
291  dest_arr(i,j,k,dest_comp) = dest_arr(i,jflip,k,dest_comp);
292  } else if (h_bc_type == ERFBCType::reflect_odd) {
293  dest_arr(i,j,k,dest_comp) = -dest_arr(i,jflip,k,dest_comp);
294  } else if (h_bc_type == ERFBCType::hoextrapcc) {
295  Real delta_j = (j - dom_hi.y);
296  dest_arr(i,j,k,dest_comp) = (1.0 + delta_j)*dest_arr(i,dom_hi.y,k,dest_comp) - delta_j*dest_arr(i,dom_hi.y-1,k,dest_comp);
297  }
298  }
299  );
300  }
301  Gpu::streamSynchronize();
302 }
RhoScalar_comp
#define RhoScalar_comp
Definition: ERF_IndexDefines.H:40
Rho_comp
#define Rho_comp
Definition: ERF_IndexDefines.H:36
RhoTheta_comp
#define RhoTheta_comp
Definition: ERF_IndexDefines.H:37
NBCVAR_max
#define NBCVAR_max
Definition: ERF_IndexDefines.H:29
NSCALARS
#define NSCALARS
Definition: ERF_IndexDefines.H:16
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)
Definition: ERF_BoundaryConditionsCons.cpp:17
BCVars::RhoScalar_bc_comp
@ RhoScalar_bc_comp
Definition: ERF_IndexDefines.H:80
ERFBCType::open
@ open
Definition: ERF_IndexDefines.H:197
ERFBCType::reflect_odd
@ reflect_odd
Definition: ERF_IndexDefines.H:187
ERFBCType::foextrap
@ foextrap
Definition: ERF_IndexDefines.H:190
ERFBCType::ext_dir
@ ext_dir
Definition: ERF_IndexDefines.H:191
ERFBCType::ext_dir_prim
@ ext_dir_prim
Definition: ERF_IndexDefines.H:193
ERFBCType::hoextrapcc
@ hoextrapcc
Definition: ERF_IndexDefines.H:198
ERFBCType::ext_dir_upwind
@ ext_dir_upwind
Definition: ERF_IndexDefines.H:199
ERFBCType::int_dir
@ int_dir
Definition: ERF_IndexDefines.H:188
ERFBCType::reflect_even
@ reflect_even
Definition: ERF_IndexDefines.H:189
MicVar_Kess::rho
@ rho
Definition: ERF_Kessler.H:22

◆ 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,
bool  do_terrain_adjustment = true 
)
321 {
322  BL_PROFILE_VAR("impose_vertical_cons_bcs()",impose_vertical_cons_bcs);
323  const auto& dom_lo = lbound(domain);
324  const auto& dom_hi = ubound(domain);
325 
326  Box per_grown_domain(domain);
327  int growx = (m_geom.isPeriodic(0)) ? 1 : 0;
328  int growy = (m_geom.isPeriodic(1)) ? 1 : 0;
329  per_grown_domain.grow(IntVect(growx,growy,0));
330  const auto& perdom_lo = lbound(per_grown_domain);
331  const auto& perdom_hi = ubound(per_grown_domain);
332 
333  GeometryData const& geomdata = m_geom.data();
334 
335  // xlo: ori = 0
336  // ylo: ori = 1
337  // zlo: ori = 2
338  // xhi: ori = 3
339  // yhi: ori = 4
340  // zhi: ori = 5
341 
342  // Based on BCRec for the domain, we need to make BCRec for this Box
343  // 0 is used as starting index for bcrs
344  Vector<BCRec> bcrs(ncomp);
345  GpuArray<GpuArray<Real, AMREX_SPACEDIM*2>,NBCVAR_max> l_bc_extdir_vals_d;
346  GpuArray<GpuArray<Real, AMREX_SPACEDIM*2>,NBCVAR_max> l_bc_neumann_vals_d;
347 
348  const int* bxlo = bx.loVect();
349  const int* bxhi = bx.hiVect();
350  const int* dlo = domain.loVect();
351  const int* dhi = domain.hiVect();
352 
353  for (int nc = 0; nc < ncomp; nc++)
354  {
355  int bc_comp = (icomp+nc >= RhoScalar_comp && icomp+nc < RhoScalar_comp+NSCALARS) ?
356  BCVars::RhoScalar_bc_comp : icomp+nc;
357  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
358  for (int dir = 0; dir < AMREX_SPACEDIM; dir++)
359  {
360  bcrs[nc].setLo(dir, ( bxlo[dir]<=dlo[dir]
361  ? m_domain_bcs_type[bc_comp].lo(dir) : BCType::int_dir ));
362  bcrs[nc].setHi(dir, ( bxhi[dir]>=dhi[dir]
363  ? m_domain_bcs_type[bc_comp].hi(dir) : BCType::int_dir ));
364  }
365 
366  for (int ori = 0; ori < 2*AMREX_SPACEDIM; ori++) {
367  l_bc_extdir_vals_d[bc_comp][ori] = m_bc_extdir_vals[bc_comp][ori];
368  l_bc_neumann_vals_d[bc_comp][ori] = m_bc_neumann_vals[bc_comp][ori];
369  }
370  }
371 
372  Gpu::DeviceVector<BCRec> bcrs_d(icomp+ncomp);
373  Gpu::copyAsync(Gpu::hostToDevice, bcrs.begin(), bcrs.end(), bcrs_d.begin());
374  const BCRec* bc_ptr = bcrs_d.data();
375 
376  {
377  Box bx_zlo(bx); bx_zlo.setBig (2,dom_lo.z-1);
378  Box bx_zhi(bx); bx_zhi.setSmall(2,dom_hi.z+1);
379  ParallelFor(
380  bx_zlo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n)
381  {
382  int dest_comp = icomp+n;
383  int bc_comp = (dest_comp >= RhoScalar_comp && dest_comp < RhoScalar_comp+NSCALARS) ?
384  BCVars::RhoScalar_bc_comp : dest_comp;
385  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
386  int l_bc_type = bc_ptr[n].lo(2);
387  if (l_bc_type == ERFBCType::ext_dir) {
388  dest_arr(i,j,k,dest_comp) = l_bc_extdir_vals_d[bc_comp][2];
389  } else if (l_bc_type == ERFBCType::ext_dir_prim) {
390  Real rho = dest_arr(i,j,dom_lo.z,Rho_comp);
391  dest_arr(i,j,k,dest_comp) = rho * l_bc_extdir_vals_d[bc_comp][2];
392  }
393  },
394  bx_zhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n)
395  {
396  int dest_comp = icomp+n;
397  int bc_comp = (dest_comp >= RhoScalar_comp && dest_comp < RhoScalar_comp+NSCALARS) ?
398  BCVars::RhoScalar_bc_comp : dest_comp;
399  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
400  int h_bc_type = bc_ptr[n].hi(2);
401  if (h_bc_type == ERFBCType::ext_dir) {
402  dest_arr(i,j,k,dest_comp) = l_bc_extdir_vals_d[bc_comp][5];
403  } else if (h_bc_type == ERFBCType::ext_dir_prim) {
404  Real rho = dest_arr(i,j,dom_hi.z,Rho_comp);
405  dest_arr(i,j,k,dest_comp) = rho * l_bc_extdir_vals_d[bc_comp][5];
406  }
407 
408  }
409  );
410  }
411 
412  {
413  Box bx_zlo(bx); bx_zlo.setBig (2,dom_lo.z-1);
414  Box bx_zhi(bx); bx_zhi.setSmall(2,dom_hi.z+1);
415  // Populate ghost cells on lo-z and hi-z domain boundaries
416  ParallelFor(
417  bx_zlo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n)
418  {
419  int dest_comp = icomp+n;
420  int bc_comp = (dest_comp >= RhoScalar_comp && dest_comp < RhoScalar_comp+NSCALARS) ?
421  BCVars::RhoScalar_bc_comp : dest_comp;
422  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
423  int l_bc_type = bc_ptr[n].lo(2);
424  int kflip = dom_lo.z - 1 - i;
425  if (l_bc_type == ERFBCType::foextrap) {
426  dest_arr(i,j,k,dest_comp) = dest_arr(i,j,dom_lo.z,dest_comp);
427  } else if (l_bc_type == ERFBCType::open) {
428  dest_arr(i,j,k,dest_comp) = dest_arr(i,j,dom_lo.z,dest_comp);
429  } else if (l_bc_type == ERFBCType::reflect_even) {
430  dest_arr(i,j,k,dest_comp) = dest_arr(i,j,kflip,dest_comp);
431  } else if (l_bc_type == ERFBCType::reflect_odd) {
432  dest_arr(i,j,k,dest_comp) = -dest_arr(i,j,kflip,dest_comp);
433  } else if (l_bc_type == ERFBCType::neumann) {
434  Real delta_z = (dom_lo.z - k) / dxInv[2];
435  dest_arr(i,j,k,dest_comp) = dest_arr(i,j,dom_lo.z,dest_comp) -
436  delta_z*l_bc_neumann_vals_d[bc_comp][2]*dest_arr(i,j,dom_lo.z,Rho_comp);
437  } else if (l_bc_type == ERFBCType::hoextrapcc) {
438  Real delta_k = (dom_lo.z - k);
439  dest_arr(i,j,k,dest_comp) = (1.0 + delta_k)*dest_arr(i,j,dom_lo.z,dest_comp) - delta_k*dest_arr(i,j,dom_lo.z+1,dest_comp);
440  }
441  },
442  bx_zhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n)
443  {
444  int dest_comp = icomp+n;
445  int bc_comp = (dest_comp >= RhoScalar_comp && dest_comp < RhoScalar_comp+NSCALARS) ?
446  BCVars::RhoScalar_bc_comp : dest_comp;
447  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
448  int h_bc_type = bc_ptr[n].hi(2);
449  int kflip = 2*dom_hi.z + 1 - i;
450  if (h_bc_type == ERFBCType::foextrap) {
451  dest_arr(i,j,k,dest_comp) = dest_arr(i,j,dom_hi.z,dest_comp);
452  } else if (h_bc_type == ERFBCType::open) {
453  dest_arr(i,j,k,dest_comp) = dest_arr(i,j,dom_hi.z,dest_comp);
454  } else if (h_bc_type == ERFBCType::reflect_even) {
455  dest_arr(i,j,k,dest_comp) = dest_arr(i,j,kflip,dest_comp);
456  } else if (h_bc_type == ERFBCType::reflect_odd) {
457  dest_arr(i,j,k,dest_comp) = -dest_arr(i,j,kflip,dest_comp);
458  } else if (h_bc_type == ERFBCType::neumann) {
459  Real delta_z = (k - dom_hi.z) / dxInv[2];
460  if( (icomp+n) == Rho_comp ) {
461  dest_arr(i,j,k,dest_comp) = dest_arr(i,j,dom_hi.z,dest_comp) +
462  delta_z*l_bc_neumann_vals_d[bc_comp][5];
463  } else {
464  dest_arr(i,j,k,dest_comp) = dest_arr(i,j,dom_hi.z,dest_comp) +
465  delta_z*l_bc_neumann_vals_d[bc_comp][5]*dest_arr(i,j,dom_hi.z,Rho_comp);
466  }
467  } else if (h_bc_type == ERFBCType::hoextrapcc){
468  Real delta_k = (k - dom_hi.z);
469  dest_arr(i,j,k,dest_comp) = (1.0 + delta_k)*dest_arr(i,j,dom_hi.z,dest_comp) - delta_k*dest_arr(i,j,dom_hi.z-1,dest_comp);
470  }
471  }
472  );
473  }
474 
475  if (do_terrain_adjustment && m_z_phys_nd) {
476  const auto& bx_lo = lbound(bx);
477  const auto& bx_hi = ubound(bx);
478  const BCRec* bc_ptr_h = bcrs.data();
479  // Neumann conditions (d<var>/dn = 0) must be aware of the surface normal with terrain.
480  // An additional source term arises from d<var>/dx & d<var>/dy & met_h_xi/eta/zeta.
481  //=====================================================================================
482  // Only modify scalars, U, or V
483  // Loop over each component
484  for (int n = 0; n < ncomp; n++) {
485  // Hit for Neumann condition at kmin
486  int dest_comp = icomp+n;
487  int l_bc_type = bc_ptr_h[n].lo(2);
488  if(l_bc_type == ERFBCType::foextrap)
489  {
490  // Loop over ghost cells in bottom XY-plane (valid box)
491  Box xybx = bx;
492 
493  int k0 = dom_lo.z;
494  if (xybx.smallEnd(2) < 0)
495  {
496  xybx.setBig(2,dom_lo.z-1);
497  xybx.setSmall(2,bx.smallEnd(2));
498 
499  // Get the dz cell size
500  Real dz = geomdata.CellSize(2);
501 
502  // Fill all the Neumann srcs with terrain
503  ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
504  {
505  // Clip indices for ghost-cells
506  int ii = amrex::min(amrex::max(i,perdom_lo.x),perdom_hi.x);
507  int jj = amrex::min(amrex::max(j,perdom_lo.y),perdom_hi.y);
508 
509  // Get metrics
510  Real met_h_xi = Compute_h_xi_AtCellCenter (ii,jj,k0,dxInv,z_phys_nd);
511  Real met_h_eta = Compute_h_eta_AtCellCenter (ii,jj,k0,dxInv,z_phys_nd);
512  Real met_h_zeta = Compute_h_zeta_AtCellCenter(ii,jj,k0,dxInv,z_phys_nd);
513 
514  // GradX at IJK location inside domain -- this relies on the assumption that we have
515  // used foextrap for cell-centered quantities outside the domain to define the gradient as zero
516  Real GradVarx, GradVary;
517  if (i < dom_lo.x-1 || i > dom_hi.x+1 || (i+1 > bx_hi.x && i-1 < bx_lo.x) ) {
518  GradVarx = 0.0;
519  } else if (i+1 > bx_hi.x) {
520  GradVarx = dxInv[0] * (dest_arr(i ,j,k0,dest_comp) - dest_arr(i-1,j,k0,dest_comp));
521  } else if (i-1 < bx_lo.x) {
522  GradVarx = dxInv[0] * (dest_arr(i+1,j,k0,dest_comp) - dest_arr(i ,j,k0,dest_comp));
523  } else {
524  GradVarx = 0.5 * dxInv[0] * (dest_arr(i+1,j,k0,dest_comp) - dest_arr(i-1,j,k0,dest_comp));
525  }
526 
527  // GradY at IJK location inside domain -- this relies on the assumption that we have
528  // used foextrap for cell-centered quantities outside the domain to define the gradient as zero
529  if (j < dom_lo.y-1 || j > dom_hi.y+1 || (j+1 > bx_hi.y && j-1 < bx_lo.y) ) {
530  GradVary = 0.0;
531  } else if (j+1 > bx_hi.y) {
532  GradVary = dxInv[1] * (dest_arr(i,j ,k0,dest_comp) - dest_arr(i,j-1,k0,dest_comp));
533  } else if (j-1 < bx_lo.y) {
534  GradVary = dxInv[1] * (dest_arr(i,j+1,k0,dest_comp) - dest_arr(i,j ,k0,dest_comp));
535  } else {
536  GradVary = 0.5 * dxInv[1] * (dest_arr(i,j+1,k0,dest_comp) - dest_arr(i,j-1,k0,dest_comp));
537  }
538 
539  // Prefactor
540  Real met_fac = met_h_zeta / ( met_h_xi*met_h_xi + met_h_eta*met_h_eta + 1. );
541 
542  // Accumulate in bottom ghost cell (EXTRAP already populated)
543  dest_arr(i,j,k,dest_comp) -= dz * met_fac * ( met_h_xi * GradVarx + met_h_eta * GradVary );
544  });
545  } // box includes k0
546  } // foextrap
547  } // ncomp
548  } // m_z_phys_nd
549  Gpu::streamSynchronize();
550 }
Compute_h_zeta_AtCellCenter
AMREX_FORCE_INLINE AMREX_GPU_DEVICE 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:46
Compute_h_eta_AtCellCenter
AMREX_GPU_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:76
Compute_h_xi_AtCellCenter
AMREX_GPU_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:61
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, bool do_terrain_adjustment=true)
Definition: ERF_BoundaryConditionsCons.cpp:317
ERFBCType::neumann
@ neumann
Definition: ERF_IndexDefines.H:195
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 amrex::Real  time,
int  bccomp_cons,
bool  do_fb = true,
bool  do_terrain_adjustment = true 
)
20 {
21  BL_PROFILE("ERFPhysBCFunct_cons::()");
22 
23  if (m_geom.isAllPeriodic()) return;
24 
25  const auto& domain = m_geom.Domain();
26  const auto dxInv = m_geom.InvCellSizeArray();
27 
28  // Create a grown domain box containing valid + periodic cells
29  Box gdomain = domain;
30  for (int i = 0; i < AMREX_SPACEDIM; ++i) {
31  if (m_geom.isPeriodic(i)) {
32  gdomain.grow(i, nghost[i]);
33  }
34  }
35 
36  //
37  // We fill all of the interior and periodic ghost cells first, so we can fill
38  // those directly inside the lateral and vertical calls.
39  //
40  if (do_fb) {
41  mf.FillBoundary(m_geom.periodicity());
42  }
43 
44 #ifdef AMREX_USE_OMP
45 #pragma omp parallel if (Gpu::notInLaunchRegion())
46 #endif
47  {
48  for (MFIter mfi(mf,false); mfi.isValid(); ++mfi)
49  {
50  //
51  // This is the box we pass to the different routines
52  // NOTE -- this is the full grid box NOT the tile box
53  //
54  Box bx = mfi.validbox();
55 
56  //
57  // These are the boxes we use to test on relative to the domain
58  //
59  Box cbx1 = bx; cbx1.grow(IntVect(nghost[0],nghost[1],0));
60  Box cbx2 = bx; cbx2.grow(nghost);
61 
62  Array4<const Real> z_nd_arr;
63 
64  if (m_z_phys_nd)
65  {
66  z_nd_arr = m_z_phys_nd->const_array(mfi);
67  }
68 
69  if (!gdomain.contains(cbx2))
70  {
71  Array4< Real> const& cons_arr = mf.array(mfi);
72  Array4<const Real> const& velx_arr = xvel.const_array(mfi);
73  Array4<const Real> const& vely_arr = yvel.const_array(mfi);
74 
75  if (!m_use_real_bcs)
76  {
77  // We send a box with ghost cells in the lateral directions only
78  impose_lateral_cons_bcs(cons_arr,velx_arr,vely_arr,cbx1,domain,icomp,ncomp,nghost);
79  }
80 
81  // We send the full FAB box with ghost cells
82  impose_vertical_cons_bcs(cons_arr,cbx2,domain,z_nd_arr,dxInv,icomp,ncomp,do_terrain_adjustment);
83  }
84 
85  } // MFIter
86  } // OpenMP
87 } // operator()
Vars::xvel
@ xvel
Definition: ERF_IndexDefines.H:141
Vars::yvel
@ yvel
Definition: ERF_IndexDefines.H:142

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_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_th_bc_data

amrex::Real* ERFPhysBCFunct_cons::m_th_bc_data
private

◆ 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: