ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
Functions
ERF_ComputeStrain_EB.cpp File Reference
#include <ERF_Diffusion.H>
#include "ERF_EddyViscosity.H"
Include dependency graph for ERF_ComputeStrain_EB.cpp:

Functions

void ComputeStrain_EB (const MFIter &mfi, Box bxcc, Box tbxxy, Box tbxxz, Box tbxyz, Box domain, const Array4< const Real > &u, const Array4< const Real > &v, const Array4< const Real > &w, Array4< Real > &tau11, Array4< Real > &tau22, Array4< Real > &tau33, Array4< Real > &tau12, Array4< Real > &tau13, Array4< Real > &tau23, const GpuArray< Real, AMREX_SPACEDIM > &dxInv, const BCRec *bc_ptr, const eb_ &ebfact, Array4< Real > &tau13i, Array4< Real > &tau23i)
 

Function Documentation

◆ ComputeStrain_EB()

void ComputeStrain_EB ( const MFIter &  mfi,
Box  bxcc,
Box  tbxxy,
Box  tbxxz,
Box  tbxyz,
Box  domain,
const Array4< const Real > &  u,
const Array4< const Real > &  v,
const Array4< const Real > &  w,
Array4< Real > &  tau11,
Array4< Real > &  tau22,
Array4< Real > &  tau33,
Array4< Real > &  tau12,
Array4< Real > &  tau13,
Array4< Real > &  tau23,
const GpuArray< Real, AMREX_SPACEDIM > &  dxInv,
const BCRec *  bc_ptr,
const eb_ &  ebfact,
Array4< Real > &  tau13i,
Array4< Real > &  tau23i 
)

Function for computing the strain rates for EB.

Parameters
[in]mfiMultiFab iterator
[in]bxcccell center box for tau_ii
[in]tbxxynodal xy box for tau_12
[in]tbxxznodal xz box for tau_13
[in]tbxyznodal yz box for tau_23
[in]domaincomputational domain
[in]ux-direction velocity
[in]vy-direction velocity
[in]wz-direction velocity
[out]tau1111 strain
[out]tau2222 strain
[out]tau3333 strain
[out]tau1212 strain
[out]tau1313 strain
[out]tau2323 strain
[in]bc_ptrcontainer with boundary condition types
[in]dxInvinverse cell size array
[in]ebfactEB factories for face-centered variables
[in]tau13icontribution to strain from du/dz
[in]tau23icontribution to strain from dv/dz
46 {
47  // Convert domain to each index type to test if we are on Dirichlet boundary
48  Box domain_xy = convert(domain, tbxxy.ixType());
49  Box domain_xz = convert(domain, tbxxz.ixType());
50  Box domain_yz = convert(domain, tbxyz.ixType());
51 
52  const auto& dom_lo = lbound(domain);
53  const auto& dom_hi = ubound(domain);
54 
55  // EB
56  // Array4<const EBCellFlag> cflag = (ebfact.get_const_factory())->getMultiEBCellFlagFab()[mfi].const_array();
57  Array4<const EBCellFlag> u_cflag = (ebfact.get_u_const_factory())->getMultiEBCellFlagFab()[mfi].const_array();
58  Array4<const EBCellFlag> v_cflag = (ebfact.get_v_const_factory())->getMultiEBCellFlagFab()[mfi].const_array();
59  Array4<const EBCellFlag> w_cflag = (ebfact.get_w_const_factory())->getMultiEBCellFlagFab()[mfi].const_array();
60 
61  // Dirichlet on left or right plane
62  bool xl_v_dir = ( (bc_ptr[BCVars::yvel_bc].lo(0) == ERFBCType::ext_dir) ||
63  (bc_ptr[BCVars::yvel_bc].lo(0) == ERFBCType::ext_dir_upwind) ||
64  (bc_ptr[BCVars::yvel_bc].lo(0) == ERFBCType::ext_dir_ingested) );
65  xl_v_dir = ( xl_v_dir && (tbxxy.smallEnd(0) == domain_xy.smallEnd(0)) );
66 
67  bool xh_v_dir = ( (bc_ptr[BCVars::yvel_bc].hi(0) == ERFBCType::ext_dir) ||
68  (bc_ptr[BCVars::yvel_bc].hi(0) == ERFBCType::ext_dir_upwind) ||
69  (bc_ptr[BCVars::yvel_bc].hi(0) == ERFBCType::ext_dir_ingested) );
70  xh_v_dir = ( xh_v_dir && (tbxxy.bigEnd(0) == domain_xy.bigEnd(0)) );
71 
72  bool xl_w_dir = ( (bc_ptr[BCVars::zvel_bc].lo(0) == ERFBCType::ext_dir) ||
73  (bc_ptr[BCVars::zvel_bc].lo(0) == ERFBCType::ext_dir_upwind) ||
74  (bc_ptr[BCVars::zvel_bc].lo(0) == ERFBCType::ext_dir_ingested) );
75  xl_w_dir = ( xl_w_dir && (tbxxz.smallEnd(0) == domain_xz.smallEnd(0)) );
76 
77  bool xh_w_dir = ( (bc_ptr[BCVars::zvel_bc].hi(0) == ERFBCType::ext_dir) ||
78  (bc_ptr[BCVars::zvel_bc].hi(0) == ERFBCType::ext_dir_upwind) ||
79  (bc_ptr[BCVars::zvel_bc].hi(0) == ERFBCType::ext_dir_ingested) );
80  xh_w_dir = ( xh_w_dir && (tbxxz.bigEnd(0) == domain_xz.bigEnd(0)) );
81 
82  // Dirichlet on front or back plane
83  bool yl_u_dir = ( (bc_ptr[BCVars::xvel_bc].lo(1) == ERFBCType::ext_dir) ||
84  (bc_ptr[BCVars::xvel_bc].lo(1) == ERFBCType::ext_dir_upwind) ||
85  (bc_ptr[BCVars::xvel_bc].lo(1) == ERFBCType::ext_dir_ingested) );
86  yl_u_dir = ( yl_u_dir && (tbxxy.smallEnd(1) == domain_xy.smallEnd(1)) );
87 
88  bool yh_u_dir = ( (bc_ptr[BCVars::xvel_bc].hi(1) == ERFBCType::ext_dir) ||
89  (bc_ptr[BCVars::xvel_bc].hi(1) == ERFBCType::ext_dir_upwind) ||
90  (bc_ptr[BCVars::xvel_bc].hi(1) == ERFBCType::ext_dir_ingested) );
91  yh_u_dir = ( yh_u_dir && (tbxxy.bigEnd(1) == domain_xy.bigEnd(1)) );
92 
93  bool yl_w_dir = ( (bc_ptr[BCVars::zvel_bc].lo(1) == ERFBCType::ext_dir) ||
94  (bc_ptr[BCVars::zvel_bc].lo(1) == ERFBCType::ext_dir_upwind) ||
95  (bc_ptr[BCVars::zvel_bc].lo(1) == ERFBCType::ext_dir_ingested) );
96  yl_w_dir = ( yl_w_dir && (tbxyz.smallEnd(1) == domain_yz.smallEnd(1)) );
97 
98  bool yh_w_dir = ( (bc_ptr[BCVars::zvel_bc].hi(1) == ERFBCType::ext_dir) ||
99  (bc_ptr[BCVars::zvel_bc].hi(1) == ERFBCType::ext_dir_upwind) ||
100  (bc_ptr[BCVars::zvel_bc].hi(1) == ERFBCType::ext_dir_ingested) );
101  yh_w_dir = ( yh_w_dir && (tbxyz.bigEnd(1) == domain_yz.bigEnd(1)) );
102 
103  // Dirichlet on top or bottom plane
104  bool zl_u_dir = ( (bc_ptr[BCVars::xvel_bc].lo(2) == ERFBCType::ext_dir) ||
105  (bc_ptr[BCVars::xvel_bc].lo(2) == ERFBCType::ext_dir_ingested) );
106  zl_u_dir = ( zl_u_dir && (tbxxz.smallEnd(2) == domain_xz.smallEnd(2)) );
107 
108  bool zh_u_dir = ( (bc_ptr[BCVars::xvel_bc].hi(2) == ERFBCType::ext_dir) ||
109  (bc_ptr[BCVars::xvel_bc].hi(2) == ERFBCType::ext_dir_ingested) );
110  zh_u_dir = ( zh_u_dir && (tbxxz.bigEnd(2) == domain_xz.bigEnd(2)) );
111 
112  bool zl_v_dir = ( (bc_ptr[BCVars::yvel_bc].lo(2) == ERFBCType::ext_dir) ||
113  (bc_ptr[BCVars::yvel_bc].lo(2) == ERFBCType::ext_dir_ingested) );
114  zl_v_dir = ( zl_v_dir && (tbxyz.smallEnd(2) == domain_yz.smallEnd(2)) );
115 
116  bool zh_v_dir = ( (bc_ptr[BCVars::yvel_bc].hi(2) == ERFBCType::ext_dir) ||
117  (bc_ptr[BCVars::yvel_bc].hi(2) == ERFBCType::ext_dir_ingested) );
118  zh_v_dir = ( zh_v_dir && (tbxyz.bigEnd(2) == domain_yz.bigEnd(2)) );
119 
120  //***********************************************************************************
121  // X-Dirichlet
122  //***********************************************************************************
123  if (xl_v_dir) {
124  Box planexy = tbxxy; planexy.setBig(0, planexy.smallEnd(0) );
125  tbxxy.growLo(0,-1);
126  bool need_to_test = (bc_ptr[BCVars::yvel_bc].lo(0) == ERFBCType::ext_dir_upwind) ? true : false;
127 
128  ParallelFor(planexy,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
129  if (!need_to_test || u(dom_lo.x,j,k) >= zero) {
130  tau12(i,j,k) = myhalf * ( (u(i, j, k) - u(i, j-1, k))*dxInv[1]
131  + (-(Real(8.)/three) * v(i-1,j,k) + three * v(i,j,k) - third * v(i+1,j,k))*dxInv[0] );
132  } else {
133  tau12(i,j,k) = myhalf * ( (u(i, j, k) - u(i, j-1, k))*dxInv[1] +
134  (v(i, j, k) - v(i-1, j, k))*dxInv[0] );
135  }
136  });
137  }
138  if (xh_v_dir) {
139  // note: tilebox xy should be nodal, so i|i-1|i-2 at the bigEnd is analogous to i-1|i|i+1 at the smallEnd
140  Box planexy = tbxxy; planexy.setSmall(0, planexy.bigEnd(0) );
141  tbxxy.growHi(0,-1);
142  bool need_to_test = (bc_ptr[BCVars::yvel_bc].hi(0) == ERFBCType::ext_dir_upwind) ? true : false;
143 
144  ParallelFor(planexy,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
145  if (!need_to_test || u(dom_hi.x+1,j,k) <= zero) {
146  tau12(i,j,k) = myhalf * ( (u(i, j, k) - u(i, j-1, k))*dxInv[1]
147  - (-(Real(8.)/three) * v(i,j,k) + three * v(i-1,j,k) - third * v(i-2,j,k))*dxInv[0] );
148  } else {
149  tau12(i,j,k) = myhalf * ( (u(i, j, k) - u(i, j-1, k))*dxInv[1] +
150  (v(i, j, k) - v(i-1, j, k))*dxInv[0] );
151  }
152  });
153  }
154 
155  if (xl_w_dir) {
156  Box planexz = tbxxz; planexz.setBig(0, planexz.smallEnd(0) );
157  tbxxz.growLo(0,-1);
158  bool need_to_test = (bc_ptr[BCVars::zvel_bc].lo(0) == ERFBCType::ext_dir_upwind) ? true : false;
159 
160  ParallelFor(planexz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
161  Real du_dz = (u(i, j, k) - u(i, j, k-1))*dxInv[2];
162  if (!need_to_test || u(dom_lo.x,j,k) >= zero) {
163  tau13(i,j,k) = myhalf * ( du_dz
164  + (-(Real(8.)/three) * w(i-1,j,k) + three * w(i,j,k) - third * w(i+1,j,k))*dxInv[0] );
165  } else {
166  tau13(i,j,k) = myhalf * ( du_dz
167  + (w(i, j, k) - w(i-1, j, k))*dxInv[0] );
168  }
169 
170  if (tau13i) tau13i(i,j,k) = myhalf * du_dz;
171  });
172  }
173  if (xh_w_dir) {
174  // note: tilebox xz should be nodal, so i|i-1|i-2 at the bigEnd is analogous to i-1|i|i+1 at the smallEnd
175  Box planexz = tbxxz; planexz.setSmall(0, planexz.bigEnd(0) );
176  tbxxz.growHi(0,-1);
177  bool need_to_test = (bc_ptr[BCVars::zvel_bc].hi(0) == ERFBCType::ext_dir_upwind) ? true : false;
178 
179  ParallelFor(planexz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
180  Real du_dz = (u(i, j, k) - u(i, j, k-1))*dxInv[2];
181  if (!need_to_test || u(dom_hi.x+1,j,k) <= zero) {
182  tau13(i,j,k) = myhalf * ( du_dz
183  - (-(Real(8.)/three) * w(i,j,k) + three * w(i-1,j,k) - third * w(i-2,j,k))*dxInv[0] );
184  } else {
185  tau13(i,j,k) = myhalf * ( du_dz
186  + (w(i, j, k) - w(i-1, j, k))*dxInv[0] );
187  }
188 
189  if (tau13i) tau13i(i,j,k) = myhalf * du_dz;
190  });
191  }
192 
193  //***********************************************************************************
194  // Y-Dirichlet
195  //***********************************************************************************
196  if (yl_u_dir) {
197  Box planexy = tbxxy; planexy.setBig(1, planexy.smallEnd(1) );
198  tbxxy.growLo(1,-1);
199  bool need_to_test = (bc_ptr[BCVars::xvel_bc].lo(1) == ERFBCType::ext_dir_upwind) ? true : false;
200 
201  ParallelFor(planexy,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
202  if (!need_to_test || v(i,dom_lo.y,k) >= zero) {
203  tau12(i,j,k) = myhalf * ( (-(Real(8.)/three) * u(i,j-1,k) + three * u(i,j,k) - third * u(i,j+1,k))*dxInv[1]
204  + (v(i, j, k) - v(i-1, j, k))*dxInv[0] );
205  } else {
206  tau12(i,j,k) = myhalf * ( (u(i, j, k) - u(i, j-1, k))*dxInv[1]
207  + (v(i, j, k) - v(i-1, j, k))*dxInv[0] );
208  }
209  });
210  }
211  if (yh_u_dir) {
212  // note: tilebox xy should be nodal, so j|j-1|j-2 at the bigEnd is analogous to j-1|j|j+1 at the smallEnd
213  Box planexy = tbxxy; planexy.setSmall(1, planexy.bigEnd(1) );
214  tbxxy.growHi(1,-1);
215  bool need_to_test = (bc_ptr[BCVars::xvel_bc].hi(1) == ERFBCType::ext_dir_upwind) ? true : false;
216 
217  ParallelFor(planexy,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
218  if (!need_to_test || v(i,dom_hi.y+1,k) <= zero) {
219  tau12(i,j,k) = myhalf * ( -(-(Real(8.)/three) * u(i,j,k) + three * u(i,j-1,k) - third * u(i,j-2,k))*dxInv[1]
220  + (v(i, j, k) - v(i-1, j, k))*dxInv[0] );
221  } else {
222  tau12(i,j,k) = myhalf * ( (u(i, j, k) - u(i, j-1, k))*dxInv[1]
223  + (v(i, j, k) - v(i-1, j, k))*dxInv[0] );
224  }
225  });
226  }
227 
228  if (yl_w_dir) {
229  Box planeyz = tbxyz; planeyz.setBig(1, planeyz.smallEnd(1) );
230  tbxyz.growLo(1,-1);
231  bool need_to_test = (bc_ptr[BCVars::zvel_bc].lo(1) == ERFBCType::ext_dir_upwind) ? true : false;
232 
233  ParallelFor(planeyz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
234  Real dv_dz = (v(i, j, k) - v(i, j, k-1))*dxInv[2];
235  if (!need_to_test || v(i,dom_lo.y,k) >= zero) {
236  tau23(i,j,k) = myhalf * ( dv_dz
237  + (-(Real(8.)/three) * w(i,j-1,k) + three * w(i,j ,k) - third * w(i,j+1,k))*dxInv[1] );
238  } else {
239  tau23(i,j,k) = myhalf * ( dv_dz
240  + (w(i, j, k) - w(i, j-1, k))*dxInv[1] );
241  }
242 
243  if (tau23i) tau23i(i,j,k) = myhalf * dv_dz;
244  });
245  }
246  if (yh_w_dir) {
247  // note: tilebox yz should be nodal, so j|j-1|j-2 at the bigEnd is analogous to j-1|j|j+1 at the smallEnd
248  Box planeyz = tbxyz; planeyz.setSmall(1, planeyz.bigEnd(1) );
249  tbxyz.growHi(1,-1);
250  bool need_to_test = (bc_ptr[BCVars::zvel_bc].hi(1) == ERFBCType::ext_dir_upwind) ? true : false;
251 
252  ParallelFor(planeyz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
253  Real dv_dz = (v(i, j, k) - v(i, j, k-1))*dxInv[2];
254  if (!need_to_test || v(i,dom_hi.y+1,k) <= zero) {
255  tau23(i,j,k) = myhalf * ( dv_dz
256  - (-(Real(8.)/three) * w(i,j ,k) + three * w(i,j-1,k) - third * w(i,j-2,k))*dxInv[1] );
257  } else {
258  tau23(i,j,k) = myhalf * ( dv_dz
259  + (w(i, j, k) - w(i, j-1, k))*dxInv[1] );
260  }
261 
262  if (tau23i) tau23i(i,j,k) = myhalf * dv_dz;
263  });
264  }
265 
266  //***********************************************************************************
267  // Z-Dirichlet
268  //***********************************************************************************
269  if (zl_u_dir) {
270  Box planexz = tbxxz; planexz.setBig(2, planexz.smallEnd(2) );
271  tbxxz.growLo(2,-1);
272 
273  ParallelFor(planexz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
274  Real du_dz = (-(Real(8.)/three) * u(i,j,k-1) + three * u(i,j,k) - third * u(i,j,k+1))*dxInv[2];
275  tau13(i,j,k) = myhalf * ( du_dz
276  + (w(i, j, k) - w(i-1, j, k))*dxInv[0] );
277 
278  if (tau13i) tau13i(i,j,k) = myhalf * du_dz;
279  });
280  }
281  if (zh_u_dir) {
282  // note: tilebox xz should be nodal, so k|k-1|k-2 at the bigEnd is analogous to k-1|k|k+1 at the smallEnd
283  Box planexz = tbxxz; planexz.setSmall(2, planexz.bigEnd(2) );
284  tbxxz.growHi(2,-1);
285 
286  ParallelFor(planexz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
287  Real du_dz = -(-(Real(8.)/three) * u(i,j,k) + three * u(i,j,k-1) - third * u(i,j,k-2))*dxInv[2];
288  tau13(i,j,k) = myhalf * ( du_dz
289  + (w(i, j, k) - w(i-1, j, k))*dxInv[0] );
290 
291  if (tau13i) tau13i(i,j,k) = myhalf * du_dz;
292  });
293  }
294 
295  if (zl_v_dir) {
296  Box planeyz = tbxyz; planeyz.setBig(2, planeyz.smallEnd(2) );
297  tbxyz.growLo(2,-1);
298 
299  ParallelFor(planeyz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
300  Real dv_dz = (-(Real(8.)/three) * v(i,j,k-1) + three * v(i,j,k ) - third * v(i,j,k+1))*dxInv[2];
301  tau23(i,j,k) = myhalf * ( dv_dz
302  + (w(i, j, k) - w(i, j-1, k))*dxInv[1] );
303 
304  if (tau23i) tau23i(i,j,k) = myhalf * dv_dz;
305  });
306  }
307  if (zh_v_dir) {
308  // note: tilebox yz should be nodal, so k|k-1|k-2 at the bigEnd is analogous to k-1|k|k+1 at the smallEnd
309  Box planeyz = tbxyz; planeyz.setSmall(2, planeyz.bigEnd(2) );
310  tbxyz.growHi(2,-1);
311 
312  ParallelFor(planeyz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
313  Real dv_dz = -(-(Real(8.)/three) * v(i,j,k ) + three * v(i,j,k-1) - third * v(i,j,k-2))*dxInv[2];
314  tau23(i,j,k) = myhalf * ( dv_dz
315  + (w(i, j, k) - w(i, j-1, k))*dxInv[1] );
316 
317  if (tau23i) tau23i(i,j,k) = myhalf * dv_dz;
318  });
319  }
320 
321  // Fill the remaining cells
322  //***********************************************************************************
323  // Cell centered strains
324  ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
325 
326  Real du_dx{zero};
327  bool can_lo_x = (i-2 >= dom_lo.x);
328  bool can_hi_x = (i+3 <= dom_hi.x+1);
329  if (can_lo_x && u_cflag(i+1,j,k).isCovered() && u_cflag(i,j,k).isSingleValued()) {
330  du_dx = ( two*u(i, j, k) - three*u(i-1, j, k) + u(i-2, j, k))*dxInv[0];
331  } else if (can_hi_x && u_cflag(i,j,k).isCovered() && u_cflag(i+1,j,k).isSingleValued()) {
332  du_dx = (- two*u(i+1, j, k) + three*u(i+2, j, k) - u(i+3, j, k))*dxInv[0];
333  } else {
334  du_dx = (u(i+1, j, k) - u(i, j, k))*dxInv[0];
335  }
336 
337  Real dv_dy{zero};
338  bool can_lo_y = (j-2 >= dom_lo.y);
339  bool can_hi_y = (j+3 <= dom_hi.y+1);
340  if (can_lo_y && v_cflag(i,j+1,k).isCovered() && v_cflag(i,j,k).isSingleValued()) {
341  dv_dy = ( two*v(i, j, k) - three*v(i, j-1, k) + v(i, j-2, k))*dxInv[1];
342  } else if (can_hi_y && v_cflag(i,j,k).isCovered() && v_cflag(i,j+1,k).isSingleValued()) {
343  dv_dy = (- two*v(i, j+1, k) + three*v(i, j+2, k) - v(i, j+3, k))*dxInv[1];
344  } else {
345  dv_dy = (v(i, j+1, k) - v(i, j, k))*dxInv[1];
346  }
347 
348  Real dw_dz{zero};
349  bool can_lo_z = (k-2 >= dom_lo.z);
350  bool can_hi_z = (k+3 <= dom_hi.z+1);
351  if (can_lo_z && w_cflag(i,j,k+1).isCovered() && w_cflag(i,j,k).isSingleValued()) {
352  dw_dz = ( two*w(i, j, k) - three*w(i, j, k-1) + w(i, j, k-2))*dxInv[2];
353  } else if (can_hi_z && w_cflag(i,j,k).isCovered() && w_cflag(i,j,k+1).isSingleValued()) {
354  dw_dz = (- two*w(i, j, k+1) + three*w(i, j, k+2) - w(i, j, k+3))*dxInv[2];
355  } else {
356  dw_dz = (w(i, j, k+1) - w(i, j, k))*dxInv[2];
357  }
358 
359  tau11(i,j,k) = du_dx;
360  tau22(i,j,k) = dv_dy;
361  tau33(i,j,k) = dw_dz;
362  });
363 
364  // Off-diagonal strains
365  ParallelFor(tbxxy,tbxxz,tbxyz,
366  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
367 
368  Real du_dy{zero};
369  bool can_lo_y_du = (j-3 >= dom_lo.y);
370  bool can_hi_y_du = (j+2 <= dom_hi.y);
371  if (can_lo_y_du && u_cflag(i,j,k).isCovered() && u_cflag(i,j-1,k).isSingleValued()) {
372  du_dy = ( two*u(i, j-1, k) - three*u(i, j-2, k) + u(i, j-3, k))*dxInv[1];
373  } else if (can_hi_y_du && u_cflag(i,j-1,k).isCovered() && u_cflag(i,j,k).isSingleValued()) {
374  du_dy = (- two*u(i, j, k) + three*u(i, j+1, k) - u(i, j+2, k))*dxInv[1];
375  } else {
376  du_dy = (u(i, j, k) - u(i, j-1, k))*dxInv[1];
377  }
378 
379  Real dv_dx{zero};
380  bool can_lo_x_dv = (i-3 >= dom_lo.x);
381  bool can_hi_x_dv = (i+2 <= dom_hi.x);
382  if (can_lo_x_dv && v_cflag(i,j,k).isCovered() && v_cflag(i-1,j,k).isSingleValued()) {
383  dv_dx = ( two*v(i-1, j, k) - three*v(i-2, j, k) + v(i-3, j, k))*dxInv[0];
384  } else if (can_hi_x_dv && v_cflag(i-1,j,k).isCovered() && v_cflag(i,j,k).isSingleValued()) {
385  dv_dx = (- two*v(i, j, k) + three*v(i+1, j, k) - v(i+2, j, k))*dxInv[0];
386  } else {
387  dv_dx = (v(i, j, k) - v(i-1, j, k))*dxInv[0];
388  }
389 
390  tau12(i,j,k) = myhalf * ( du_dy + dv_dx );
391  },
392  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
393 
394  Real du_dz{zero};
395  bool can_lo_z_du = (k-3 >= dom_lo.z);
396  bool can_hi_z_du = (k+2 <= dom_hi.z);
397  if (can_lo_z_du && u_cflag(i,j,k).isCovered() && u_cflag(i,j,k-1).isSingleValued()) {
398  du_dz = ( two*u(i, j, k-1) - three*u(i, j, k-2) + u(i, j, k-3))*dxInv[2];
399  } else if (can_hi_z_du && u_cflag(i,j,k-1).isCovered() && u_cflag(i,j,k).isSingleValued()) {
400  du_dz = (- two*u(i, j, k) + three*u(i, j, k+1) - u(i, j, k+2))*dxInv[2];
401  } else {
402  du_dz = (u(i, j, k) - u(i, j, k-1))*dxInv[2];
403  }
404 
405  Real dw_dx{zero};
406  bool can_lo_x_dw = (i-3 >= dom_lo.x);
407  bool can_hi_x_dw = (i+2 <= dom_hi.x);
408  if (can_lo_x_dw && w_cflag(i,j,k).isCovered() && w_cflag(i-1,j,k).isSingleValued()) {
409  dw_dx = ( two*w(i-1, j, k) - three*w(i-2, j, k) + w(i-3, j, k))*dxInv[0];
410  } else if (can_hi_x_dw && w_cflag(i-1,j,k).isCovered() && w_cflag(i,j,k).isSingleValued()) {
411  dw_dx = (- two*w(i, j, k) + three*w(i+1, j, k) - w(i+2, j, k))*dxInv[0];
412  } else {
413  dw_dx = (w(i, j, k) - w(i-1, j, k))*dxInv[0];
414  }
415 
416  tau13(i,j,k) = myhalf * ( du_dz + dw_dx );
417 
418  if (tau13i) tau13i(i,j,k) = myhalf * du_dz;
419  },
420  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
421 
422  Real dv_dz{zero};
423  bool can_lo_z_dv = (k-3 >= dom_lo.z);
424  bool can_hi_z_dv = (k+2 <= dom_hi.z);
425  if (can_lo_z_dv && v_cflag(i,j,k).isCovered() && v_cflag(i,j,k-1).isSingleValued()) {
426  dv_dz = ( two*v(i, j, k-1) - three*v(i, j, k-2) + v(i, j, k-3))*dxInv[2];
427  } else if (can_hi_z_dv && v_cflag(i,j,k-1).isCovered() && v_cflag(i,j,k).isSingleValued()) {
428  dv_dz = (- two*v(i, j, k) + three*v(i, j, k+1) - v(i, j, k+2))*dxInv[2];
429  } else {
430  dv_dz = (v(i, j, k) - v(i, j, k-1))*dxInv[2];
431  }
432 
433  Real dw_dy{zero};
434  bool can_lo_y_dw = (j-3 >= dom_lo.y);
435  bool can_hi_y_dw = (j+2 <= dom_hi.y);
436  if (can_lo_y_dw && w_cflag(i,j,k).isCovered() && w_cflag(i,j-1,k).isSingleValued()) {
437  dw_dy = ( two*w(i, j-1, k) - three*w(i, j-2, k) + w(i, j-3, k))*dxInv[1];
438  } else if (can_hi_y_dw && w_cflag(i,j-1,k).isCovered() && w_cflag(i,j,k).isSingleValued()) {
439  dw_dy = (- two*w(i, j, k) + three*w(i, j+1, k) - w(i, j+2, k))*dxInv[1];
440  } else {
441  dw_dy = (w(i, j, k) - w(i, j-1, k))*dxInv[1];
442  }
443 
444 
445  tau23(i,j,k) = myhalf * ( dv_dz + dw_dy );
446 
447  if (tau23i) tau23i(i,j,k) = myhalf * dv_dz;
448  });
449 }
three
constexpr amrex::Real three
Definition: ERF_Constants.H:11
two
constexpr amrex::Real two
Definition: ERF_Constants.H:10
zero
constexpr amrex::Real zero
Definition: ERF_Constants.H:8
myhalf
constexpr amrex::Real myhalf
Definition: ERF_Constants.H:13
third
constexpr amrex::Real third
Definition: ERF_Constants.H:15
tau12
@ tau12
Definition: ERF_DataStruct.H:38
tau23
@ tau23
Definition: ERF_DataStruct.H:38
tau33
@ tau33
Definition: ERF_DataStruct.H:38
tau22
@ tau22
Definition: ERF_DataStruct.H:38
tau11
@ tau11
Definition: ERF_DataStruct.H:38
tau13
@ tau13
Definition: ERF_DataStruct.H:38
dxInv
amrex::GpuArray< Real, AMREX_SPACEDIM > dxInv
Definition: ERF_InitCustomPertVels_ParticleTests.H:17
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);})
w
Real w
Definition: ERF_Plotfile2DInterpolator.cpp:22
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:19
eb_::get_w_const_factory
eb_aux_ const * get_w_const_factory() const noexcept
Return the ERF auxiliary z-face EB factory.
Definition: ERF_EB.H:123
eb_::get_v_const_factory
eb_aux_ const * get_v_const_factory() const noexcept
Return the ERF auxiliary y-face EB factory.
Definition: ERF_EB.H:121
eb_::get_u_const_factory
eb_aux_ const * get_u_const_factory() const noexcept
Return the ERF auxiliary x-face EB factory.
Definition: ERF_EB.H:119
BCVars::zvel_bc
@ zvel_bc
Definition: ERF_IndexDefines.H:104
BCVars::yvel_bc
@ yvel_bc
Definition: ERF_IndexDefines.H:103
BCVars::xvel_bc
@ xvel_bc
Definition: ERF_IndexDefines.H:102
ERFBCType::ext_dir_ingested
@ ext_dir_ingested
Definition: ERF_IndexDefines.H:253
ERFBCType::ext_dir
@ ext_dir
Definition: ERF_IndexDefines.H:249
ERFBCType::ext_dir_upwind
@ ext_dir_upwind
Definition: ERF_IndexDefines.H:257

Referenced by ERF::advance_dycore(), and erf_make_tau_terms().

Here is the call graph for this function:
Here is the caller graph for this function: