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

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: