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

Functions

void ComputeStrain_T (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 > &tau21, Array4< Real > &tau13, Array4< Real > &tau31, Array4< Real > &tau23, Array4< Real > &tau32, const Array4< const Real > &z_nd, const Array4< const Real > &detJ, const GpuArray< Real, AMREX_SPACEDIM > &dxInv, const Array4< const Real > &mf_mx, const Array4< const Real > &mf_ux, const Array4< const Real > &mf_vx, const Array4< const Real > &mf_my, const Array4< const Real > &mf_uy, const Array4< const Real > &mf_vy, const BCRec *bc_ptr, Array4< Real > &tau13i, Array4< Real > &tau23i)
 

Function Documentation

◆ ComputeStrain_T()

void ComputeStrain_T ( 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 > &  tau21,
Array4< Real > &  tau13,
Array4< Real > &  tau31,
Array4< Real > &  tau23,
Array4< Real > &  tau32,
const Array4< const Real > &  z_nd,
const Array4< const Real > &  detJ,
const GpuArray< Real, AMREX_SPACEDIM > &  dxInv,
const Array4< const Real > &  mf_mx,
const Array4< const Real > &  mf_ux,
const Array4< const Real > &  mf_vx,
const Array4< const Real > &  mf_my,
const Array4< const Real > &  mf_uy,
const Array4< const Real > &  mf_vy,
const BCRec *  bc_ptr,
Array4< Real > &  tau13i,
Array4< Real > &  tau23i 
)

Function for computing the strain rates with terrain.

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]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]tau2121 strain
[out]tau2323 strain
[out]tau3131 strain
[out]tau3232 strain
[in]z_ndnodal array of physical z heights
[in]detJJacobian determinant
[in]bc_ptrcontainer with boundary condition types
[in]dxInvinverse cell size array
[in]mf_mxmap factor at cell center
[in]mf_uxmap factor at x-face
[in]mf_vxmap factor at y-face
[in]mf_mymap factor at cell center
[in]mf_uymap factor at x-face
[in]mf_vymap factor at y-face
[in]tau13icontribution to strain from du/dz
[in]tau23icontribution to strain from dv/dz
62 {
63  // Convert domain to each index type to test if we are on dirichlet boundary
64  Box domain_xy = convert(domain, tbxxy.ixType());
65  Box domain_xz = convert(domain, tbxxz.ixType());
66  Box domain_yz = convert(domain, tbxyz.ixType());
67 
68  const auto& dom_lo = lbound(domain);
69  const auto& dom_hi = ubound(domain);
70 
71  // Dirichlet on left or right plane
72  bool xl_v_dir = ( (bc_ptr[BCVars::yvel_bc].lo(0) == ERFBCType::ext_dir) ||
73  (bc_ptr[BCVars::yvel_bc].lo(0) == ERFBCType::ext_dir_upwind) ||
74  (bc_ptr[BCVars::yvel_bc].lo(0) == ERFBCType::ext_dir_ingested) );
75  xl_v_dir = ( xl_v_dir && (tbxxy.smallEnd(0) == domain_xy.smallEnd(0)) );
76 
77  bool xh_v_dir = ( (bc_ptr[BCVars::yvel_bc].hi(0) == ERFBCType::ext_dir) ||
78  (bc_ptr[BCVars::yvel_bc].hi(0) == ERFBCType::ext_dir_upwind) ||
79  (bc_ptr[BCVars::yvel_bc].hi(0) == ERFBCType::ext_dir_ingested) );
80  xh_v_dir = ( xh_v_dir && (tbxxy.bigEnd(0) == domain_xy.bigEnd(0)) );
81 
82  bool xl_w_dir = ( (bc_ptr[BCVars::zvel_bc].lo(0) == ERFBCType::ext_dir) ||
83  (bc_ptr[BCVars::zvel_bc].lo(0) == ERFBCType::ext_dir_upwind) ||
84  (bc_ptr[BCVars::zvel_bc].lo(0) == ERFBCType::ext_dir_ingested) );
85  xl_w_dir = ( xl_w_dir && (tbxxz.smallEnd(0) == domain_xz.smallEnd(0)) );
86 
87  bool xh_w_dir = ( (bc_ptr[BCVars::zvel_bc].hi(0) == ERFBCType::ext_dir) ||
88  (bc_ptr[BCVars::zvel_bc].hi(0) == ERFBCType::ext_dir_upwind) ||
89  (bc_ptr[BCVars::zvel_bc].hi(0) == ERFBCType::ext_dir_ingested) );
90  xh_w_dir = ( xh_w_dir && (tbxxz.bigEnd(0) == domain_xz.bigEnd(0)) );
91 
92  // Dirichlet on front or back plane
93  bool yl_u_dir = ( (bc_ptr[BCVars::xvel_bc].lo(1) == ERFBCType::ext_dir) ||
94  (bc_ptr[BCVars::xvel_bc].lo(1) == ERFBCType::ext_dir_upwind) ||
95  (bc_ptr[BCVars::xvel_bc].lo(1) == ERFBCType::ext_dir_ingested) );
96  yl_u_dir = ( yl_u_dir && (tbxxy.smallEnd(1) == domain_xy.smallEnd(1)) );
97 
98  bool yh_u_dir = ( (bc_ptr[BCVars::xvel_bc].hi(1) == ERFBCType::ext_dir) ||
99  (bc_ptr[BCVars::xvel_bc].hi(1) == ERFBCType::ext_dir_upwind) ||
100  (bc_ptr[BCVars::xvel_bc].hi(1) == ERFBCType::ext_dir_ingested) );
101  yh_u_dir = ( yh_u_dir && (tbxxy.bigEnd(1) == domain_xy.bigEnd(1)) );
102 
103  bool yl_w_dir = ( (bc_ptr[BCVars::zvel_bc].lo(1) == ERFBCType::ext_dir) ||
104  (bc_ptr[BCVars::zvel_bc].lo(1) == ERFBCType::ext_dir_upwind) ||
105  (bc_ptr[BCVars::zvel_bc].lo(1) == ERFBCType::ext_dir_ingested) );
106  yl_w_dir = ( yl_w_dir && (tbxyz.smallEnd(1) == domain_yz.smallEnd(1)) );
107 
108  bool yh_w_dir = ( (bc_ptr[BCVars::zvel_bc].hi(1) == ERFBCType::ext_dir) ||
109  (bc_ptr[BCVars::zvel_bc].hi(1) == ERFBCType::ext_dir_upwind) ||
110  (bc_ptr[BCVars::zvel_bc].hi(1) == ERFBCType::ext_dir_ingested) );
111  yh_w_dir = ( yh_w_dir && (tbxyz.bigEnd(1) == domain_yz.bigEnd(1)) );
112 
113  // Dirichlet on top or bottom plane
114  bool zl_u_dir = ( (bc_ptr[BCVars::xvel_bc].lo(2) == ERFBCType::ext_dir) ||
115  (bc_ptr[BCVars::xvel_bc].lo(2) == ERFBCType::ext_dir_ingested) );
116  zl_u_dir = ( zl_u_dir && (tbxxz.smallEnd(2) == domain_xz.smallEnd(2)) );
117 
118  bool zh_u_dir = ( (bc_ptr[BCVars::xvel_bc].hi(2) == ERFBCType::ext_dir) ||
119  (bc_ptr[BCVars::xvel_bc].hi(2) == ERFBCType::ext_dir_ingested) );
120  zh_u_dir = ( zh_u_dir && (tbxxz.bigEnd(2) == domain_xz.bigEnd(2)) );
121 
122  bool zl_v_dir = ( (bc_ptr[BCVars::yvel_bc].lo(2) == ERFBCType::ext_dir) ||
123  (bc_ptr[BCVars::yvel_bc].lo(2) == ERFBCType::ext_dir_ingested) );
124  zl_v_dir = ( zl_v_dir && (tbxyz.smallEnd(2) == domain_yz.smallEnd(2)) );
125 
126  bool zh_v_dir = ( (bc_ptr[BCVars::yvel_bc].hi(2) == ERFBCType::ext_dir) ||
127  (bc_ptr[BCVars::yvel_bc].hi(2) == ERFBCType::ext_dir_ingested) );
128  zh_v_dir = ( zh_v_dir && (tbxyz.bigEnd(2) == domain_yz.bigEnd(2)) );
129 
130  //***********************************************************************************
131  // X-Dirichlet
132  //***********************************************************************************
133  if (xl_v_dir) {
134  Box planexy = tbxxy; planexy.setBig(0, planexy.smallEnd(0) );
135  tbxxy.growLo(0,-1);
136  bool need_to_test = (bc_ptr[BCVars::yvel_bc].lo(0) == ERFBCType::ext_dir_upwind) ? true : false;
137 
138  ParallelFor(planexy,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
139  Real inv_dist = (k == 0) ? two / (z_nd(i,j,k+2) - z_nd(i,j,k )) :
140  two / (z_nd(i,j,k+2) + z_nd(i,j,k+1) - z_nd(i,j,k) - z_nd(i,j,k-1));
141 
142  Real GradUz = (k == 0) ?
143  inv_dist * ( u(i ,j ,k+1) + u(i ,j-1,k+1)
144  - u(i ,j ,k ) - u(i ,j-1,k ) ) :
145  inv_dist * ( u(i ,j ,k+1) + u(i ,j-1,k+1)
146  - u(i ,j ,k-1) - u(i ,j-1,k-1) );
147  Real GradVz = (k == 0) ?
148  inv_dist * ( v(i ,j ,k+1) + v(i-1,j ,k+1)
149  - v(i ,j ,k ) - v(i-1,j ,k ) ) :
150  inv_dist * ( v(i ,j ,k+1) + v(i-1,j ,k+1)
151  - v(i ,j ,k-1) - v(i-1,j ,k-1) );
152 
153  Real mfy = myhalf * (mf_uy(i,j,0) + mf_uy(i ,j-1,0));
154  Real mfx = myhalf * (mf_vx(i,j,0) + mf_vx(i-1,j ,0));
155 
156  Real met_h_xi = Compute_h_xi_AtEdgeCenterK (i,j,k,dxInv,z_nd);
157  Real met_h_eta = Compute_h_eta_AtEdgeCenterK (i,j,k,dxInv,z_nd);
158 
159  if (!need_to_test || u(dom_lo.x,j,k) >= zero) {
160  tau12(i,j,k) = myhalf * ( (u(i, j, k) - u(i, j-1, k) )*dxInv[1]*mfy
161  + (-(Real(8.)/three) * v(i-1,j,k) + three * v(i,j,k) - (one/three) * v(i+1,j,k))*dxInv[0]*mfx
162  - (met_h_eta)*GradUz*mfy
163  - (met_h_xi )*GradVz*mfx );
164  // Note:
165  // tau12 ~ (du/dη) * (dη/dy)
166  // + (dv/dξ) * (dξ/dx)
167  // - (dz/dη) * (du/dz) * (dη/dy)
168  // - (dz/dξ) * (dv/dz) * (dξ/dx)
169  // ~ du/dy + dv/dx
170  // where ξ and η are in computational space, corresponding to
171  // x and y in physical space
172  } else {
173  tau12(i,j,k) = myhalf * ( (u(i, j, k) - u(i , j-1, k))*dxInv[1]*mfy
174  + (v(i, j, k) - v(i-1, j , k))*dxInv[0]*mfx
175  - (met_h_eta)*GradUz*mfy
176  - (met_h_xi )*GradVz*mfx );
177  }
178  tau21(i,j,k) = tau12(i,j,k);
179  });
180  }
181  if (xh_v_dir) {
182  Box planexy = tbxxy; planexy.setSmall(0, planexy.bigEnd(0) );
183  tbxxy.growHi(0,-1);
184  bool need_to_test = (bc_ptr[BCVars::yvel_bc].hi(0) == ERFBCType::ext_dir_upwind) ? true : false;
185 
186  ParallelFor(planexy,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
187  Real inv_dist = (k == 0) ? two / (z_nd(i,j,k+2) - z_nd(i,j,k )) :
188  two / (z_nd(i,j,k+2) + z_nd(i,j,k+1) - z_nd(i,j,k) - z_nd(i,j,k-1));
189 
190  Real GradUz = (k == 0) ?
191  inv_dist * ( u(i ,j ,k+1) + u(i ,j-1,k+1)
192  - u(i ,j ,k ) - u(i ,j-1,k ) ) :
193  inv_dist * ( u(i ,j ,k+1) + u(i ,j-1,k+1)
194  - u(i ,j ,k-1) - u(i ,j-1,k-1) );
195  Real GradVz = (k == 0) ?
196  inv_dist * ( v(i ,j ,k+1) + v(i-1,j ,k+1)
197  - v(i ,j ,k ) - v(i-1,j ,k ) ) :
198  inv_dist * ( v(i ,j ,k+1) + v(i-1,j ,k+1)
199  - v(i ,j ,k-1) - v(i-1,j ,k-1) );
200 
201  Real mfy = myhalf * (mf_uy(i,j,0) + mf_uy(i ,j-1,0));
202  Real mfx = myhalf * (mf_vx(i,j,0) + mf_vx(i-1,j ,0));
203 
204  Real met_h_xi = Compute_h_xi_AtEdgeCenterK (i,j,k,dxInv,z_nd);
205  Real met_h_eta = Compute_h_eta_AtEdgeCenterK (i,j,k,dxInv,z_nd);
206 
207  if (!need_to_test || u(dom_hi.x+1,j,k) <= zero) {
208  tau12(i,j,k) = myhalf * ( (u(i, j, k) - u(i, j-1, k) )*dxInv[1]*mfy
209  - (-(Real(8.)/three) * v(i,j,k) + three * v(i-1,j,k) - (one/three) * v(i-2,j,k))*dxInv[0]*mfx
210  - (met_h_eta)*GradUz*mfy
211  - (met_h_xi )*GradVz*mfx );
212  } else {
213  tau12(i,j,k) = myhalf * ( (u(i, j, k) - u(i , j-1, k))*dxInv[1]*mfy
214  + (v(i, j, k) - v(i-1, j , k))*dxInv[0]*mfx
215  - (met_h_eta)*GradUz*mfy
216  - (met_h_xi )*GradVz*mfx );
217  }
218  tau21(i,j,k) = tau12(i,j,k);
219  });
220  }
221 
222  if (xl_w_dir) {
223  Box planexz = tbxxz; planexz.setBig(0, planexz.smallEnd(0) );
224  planexz.setSmall(2, planexz.smallEnd(2)+1 ); planexz.setBig(2, planexz.bigEnd(2)-1 );
225  tbxxz.growLo(0,-1);
226  bool need_to_test = (bc_ptr[BCVars::zvel_bc].lo(0) == ERFBCType::ext_dir_upwind) ? true : false;
227 
228  ParallelFor(planexz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
229  Real dz0 = myhalf * ( z_nd(i,j,k+1) + z_nd(i,j+1,k+1)
230  - z_nd(i,j,k-1) - z_nd(i,j+1,k-1) );
231  Real idz0 = one / dz0;
232 
233  Real GradWz = myhalf * idz0 * ( w(i ,j ,k+1) + w(i-1,j ,k+1)
234  - w(i ,j ,k-1) - w(i-1,j ,k-1) );
235  Real mfx = mf_ux(i,j,0);
236 
237  Real met_h_xi = Compute_h_xi_AtEdgeCenterJ (i,j,k,dxInv,z_nd);
238  Real met_h_zeta = Compute_h_zeta_AtEdgeCenterJ(i,j,k,dxInv,z_nd);
239 
240  Real du_dz = (u(i, j, k) - u(i, j, k-1))*dxInv[2]/met_h_zeta;
241  if (!need_to_test || u(dom_lo.x,j,k) >= zero) {
242  tau13(i,j,k) = myhalf * ( du_dz
243  + ( (-(Real(8.)/three) * w(i-1,j,k) + three * w(i,j,k) - (one/three) * w(i+1,j,k))*dxInv[0]
244  - (met_h_xi)*GradWz ) * mfx );
245  // Note:
246  // tau13 ~ (du/dζ) / (dz/dζ)
247  // + (dw/dξ) * (dξ/dx)
248  // - (dz/dξ) * (dw/dz) * (dξ/dx)
249  // ~ du/dz + dw/dx
250  // where ξ and ζ are in computational space, corresponding to
251  // x and z in physical space
252  } else {
253  tau13(i,j,k) = myhalf * ( du_dz
254  + ( (w(i, j, k) - w(i-1, j, k ))*dxInv[0]
255  - (met_h_xi)*GradWz ) * mfx );
256  }
257  tau31(i,j,k) = tau13(i,j,k);
258 
259  if (tau13i) tau13i(i,j,k) = myhalf * du_dz;
260  });
261  }
262 
263  if (xh_w_dir) {
264  Box planexz = tbxxz; planexz.setSmall(0, planexz.bigEnd(0) );
265  planexz.setSmall(2, planexz.smallEnd(2)+1 ); planexz.setBig(2, planexz.bigEnd(2)-1 );
266  tbxxz.growHi(0,-1);
267  bool need_to_test = (bc_ptr[BCVars::zvel_bc].hi(0) == ERFBCType::ext_dir_upwind) ? true : false;
268 
269  ParallelFor(planexz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
270  Real dz0 = myhalf * ( z_nd(i,j,k+1) + z_nd(i,j+1,k+1)
271  - z_nd(i,j,k-1) - z_nd(i,j+1,k-1) );
272  Real idz0 = one / dz0;
273 
274  Real GradWz = myhalf * idz0 * ( w(i ,j ,k+1) + w(i-1,j ,k+1)
275  - w(i ,j ,k-1) - w(i-1,j ,k-1) );
276 
277  Real mfx = mf_ux(i,j,0);
278 
279  Real met_h_xi = Compute_h_xi_AtEdgeCenterJ (i,j,k,dxInv,z_nd);
280  Real met_h_zeta = Compute_h_zeta_AtEdgeCenterJ(i,j,k,dxInv,z_nd);
281 
282  Real du_dz = (u(i, j, k) - u(i, j, k-1))*dxInv[2]/met_h_zeta;
283  if (!need_to_test || u(dom_hi.x+1,j,k) <= zero) {
284  tau13(i,j,k) = myhalf * ( du_dz
285  - ( (-(Real(8.)/three) * w(i,j,k) + three * w(i-1,j,k) - (one/three) * w(i-2,j,k))*dxInv[0]
286  - (met_h_xi)*GradWz ) * mfx );
287  } else {
288  tau13(i,j,k) = myhalf * ( du_dz
289  + ( (w(i, j, k) - w(i-1, j, k ))*dxInv[0]
290  - (met_h_xi)*GradWz ) * mfx );
291  }
292  tau31(i,j,k) = tau13(i,j,k);
293 
294  if (tau13i) tau13i(i,j,k) = myhalf * du_dz;
295  });
296  }
297 
298  //***********************************************************************************
299  // Y-Dirichlet
300  //***********************************************************************************
301  if (yl_u_dir) {
302  Box planexy = tbxxy; planexy.setBig(1, planexy.smallEnd(1) );
303  tbxxy.growLo(1,-1);
304  bool need_to_test = (bc_ptr[BCVars::xvel_bc].lo(1) == ERFBCType::ext_dir_upwind) ? true : false;
305 
306  ParallelFor(planexy,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
307  Real dz0 = ( z_nd(i,j,k+1) - z_nd(i,j,k-1) );
308  Real idz0 = one / dz0;
309 
310  Real GradUz = (k == 0) ?
311  idz0 * ( u(i ,j ,k+1) + u(i ,j-1,k+1)
312  - u(i ,j ,k ) - u(i ,j-1,k ) ) :
313  myhalf * idz0 * ( u(i ,j ,k+1) + u(i ,j-1,k+1)
314  - u(i ,j ,k-1) - u(i ,j-1,k-1) );
315  Real GradVz = (k == 0) ?
316  idz0 * ( v(i ,j ,k+1) + v(i-1,j ,k+1)
317  - v(i ,j ,k ) - v(i-1,j ,k ) ) :
318  myhalf * idz0 * ( v(i ,j ,k+1) + v(i-1,j ,k+1)
319  - v(i ,j ,k-1) - v(i-1,j ,k-1) );
320 
321  Real mfy = myhalf * (mf_uy(i,j,0) + mf_uy(i ,j-1,0));
322  Real mfx = myhalf * (mf_vx(i,j,0) + mf_vx(i-1,j ,0));
323 
324  Real met_h_xi = Compute_h_xi_AtEdgeCenterK (i,j,k,dxInv,z_nd);
325  Real met_h_eta = Compute_h_eta_AtEdgeCenterK (i,j,k,dxInv,z_nd);
326 
327  if (!need_to_test || v(i,dom_lo.y,k) >= zero) {
328  tau12(i,j,k) = myhalf * ( (-(Real(8.)/three) * u(i,j-1,k) + three * u(i,j,k) - (one/three) * u(i,j+1,k))*dxInv[1]*mfy
329  + (v(i, j, k) - v(i-1, j, k))*dxInv[0]*mfx
330  - (met_h_eta)*GradUz*mfy
331  - (met_h_xi )*GradVz*mfx );
332  } else {
333  tau12(i,j,k) = myhalf * ( (u(i, j, k) - u(i , j-1, k))*dxInv[1]*mfy
334  + (v(i, j, k) - v(i-1, j , k))*dxInv[0]*mfx
335  - (met_h_eta)*GradUz*mfy
336  - (met_h_xi )*GradVz*mfx );
337  }
338  tau21(i,j,k) = tau12(i,j,k);
339  });
340  }
341  if (yh_u_dir) {
342  Box planexy = tbxxy; planexy.setSmall(1, planexy.bigEnd(1) );
343  tbxxy.growHi(1,-1);
344  bool need_to_test = (bc_ptr[BCVars::xvel_bc].hi(1) == ERFBCType::ext_dir_upwind) ? true : false;
345 
346  ParallelFor(planexy,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
347  Real dz0 = ( z_nd(i,j,k+1) - z_nd(i,j,k-1) );
348  Real idz0 = one / dz0;
349 
350  Real GradUz = (k == 0) ?
351  idz0 * ( u(i ,j ,k+1) + u(i ,j-1,k+1)
352  - u(i ,j ,k ) - u(i ,j-1,k ) ) :
353  myhalf * idz0 * ( u(i ,j ,k+1) + u(i ,j-1,k+1)
354  - u(i ,j ,k-1) - u(i ,j-1,k-1) );
355  Real GradVz = (k == 0) ?
356  idz0 * ( v(i ,j ,k+1) + v(i-1,j ,k+1)
357  - v(i ,j ,k ) - v(i-1,j ,k ) ) :
358  myhalf * idz0 * ( v(i ,j ,k+1) + v(i-1,j ,k+1)
359  - v(i ,j ,k-1) - v(i-1,j ,k-1) );
360 
361  Real mfy = myhalf * (mf_uy(i,j,0) + mf_uy(i ,j-1,0));
362  Real mfx = myhalf * (mf_vx(i,j,0) + mf_vx(i-1,j ,0));
363 
364  Real met_h_xi = Compute_h_xi_AtEdgeCenterK (i,j,k,dxInv,z_nd);
365  Real met_h_eta = Compute_h_eta_AtEdgeCenterK (i,j,k,dxInv,z_nd);
366 
367  if (!need_to_test || v(i,dom_hi.y+1,k) <= zero) {
368  tau12(i,j,k) = myhalf * ( -(-(Real(8.)/three) * u(i,j,k) + three * u(i,j-1,k) - (one/three) * u(i,j-2,k))*dxInv[1]*mfy +
369  + (v(i, j, k) - v(i-1, j, k))*dxInv[0]*mfx
370  - (met_h_eta)*GradUz*mfy
371  - (met_h_xi )*GradVz*mfx );
372  } else {
373  tau12(i,j,k) = myhalf * ( (u(i, j, k) - u(i , j-1, k))*dxInv[1]*mfy
374  + (v(i, j, k) - v(i-1, j , k))*dxInv[0]*mfx
375  - (met_h_eta)*GradUz*mfy
376  - (met_h_xi )*GradVz*mfx );
377  }
378  tau21(i,j,k) = tau12(i,j,k);
379  });
380  }
381 
382  if (yl_w_dir) {
383  Box planeyz = tbxyz; planeyz.setBig(1, planeyz.smallEnd(1) );
384  planeyz.setSmall(2, planeyz.smallEnd(2)+1 ); planeyz.setBig(2, planeyz.bigEnd(2)-1 );
385  tbxyz.growLo(1,-1);
386  bool need_to_test = (bc_ptr[BCVars::zvel_bc].lo(1) == ERFBCType::ext_dir_upwind) ? true : false;
387 
388  ParallelFor(planeyz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
389  Real dz0 = myhalf * ( z_nd(i,j,k+1) + z_nd(i+1,j,k+1)
390  - z_nd(i,j,k-1) - z_nd(i+1,j,k-1) );
391  Real idz0 = one / dz0;
392 
393  Real GradWz = myhalf * idz0 * ( w(i ,j ,k+1) + w(i ,j-1,k+1)
394  - w(i ,j ,k-1) - w(i ,j-1,k-1) );
395 
396  Real mfy = mf_vy(i,j,0);
397 
398  Real met_h_eta = Compute_h_eta_AtEdgeCenterI (i,j,k,dxInv,z_nd);
399  Real met_h_zeta = Compute_h_zeta_AtEdgeCenterI(i,j,k,dxInv,z_nd);
400 
401  Real dv_dz = (v(i, j, k) - v(i, j, k-1))*dxInv[2]/met_h_zeta;
402  if (!need_to_test || v(i,dom_lo.y,k) >= zero) {
403  tau23(i,j,k) = myhalf * ( dv_dz
404  + ( (-(Real(8.)/three) * w(i,j-1,k) + three * w(i,j ,k) - (one/three) * w(i,j+1,k))*dxInv[1]
405  - (met_h_eta)*GradWz ) * mfy );
406  // Note:
407  // tau23 ~ (dv/dζ) / (dz/dζ)
408  // + (dw/dη) * (dη/dy)
409  // - (dz/dη) * (dw/dz) * (dη/dy)
410  // ~ dv/dz + dw/dy
411  // where η and ζ are in computational space, corresponding to
412  // y and z in physical space
413  } else {
414  tau23(i,j,k) = myhalf * ( dv_dz
415  + ( (w(i, j, k) - w(i, j-1, k ))*dxInv[1]
416  - (met_h_eta)*GradWz ) * mfy );
417  }
418  tau32(i,j,k) = tau23(i,j,k);
419 
420  if (tau23i) tau23i(i,j,k) = myhalf * dv_dz;
421  });
422  }
423  if (yh_w_dir) {
424  Box planeyz = tbxyz; planeyz.setSmall(1, planeyz.bigEnd(1) );
425  planeyz.setSmall(2, planeyz.smallEnd(2)+1 ); planeyz.setBig(2, planeyz.bigEnd(2)-1 );
426  tbxyz.growHi(1,-1);
427  bool need_to_test = (bc_ptr[BCVars::zvel_bc].hi(1) == ERFBCType::ext_dir_upwind) ? true : false;
428 
429  ParallelFor(planeyz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
430  Real dz0 = myhalf * ( z_nd(i,j,k+1) + z_nd(i+1,j,k+1)
431  - z_nd(i,j,k-1) - z_nd(i+1,j,k-1) );
432  Real idz0 = one / dz0;
433 
434  Real GradWz = myhalf * idz0 * ( w(i ,j ,k+1) + w(i ,j-1,k+1)
435  - w(i ,j ,k-1) - w(i ,j-1,k-1) );
436 
437  Real mfy = mf_vy(i,j,0);
438 
439  Real met_h_eta = Compute_h_eta_AtEdgeCenterI (i,j,k,dxInv,z_nd);
440  Real met_h_zeta = Compute_h_zeta_AtEdgeCenterI(i,j,k,dxInv,z_nd);
441 
442  Real dv_dz = (v(i, j, k) - v(i, j, k-1))*dxInv[2]/met_h_zeta;
443  if (!need_to_test || v(i,dom_hi.y+1,k) <= zero) {
444  tau23(i,j,k) = myhalf * ( dv_dz
445  - ( (-(Real(8.)/three) * w(i,j ,k) + three * w(i,j-1,k) - (one/three) * w(i,j-2,k))*dxInv[1]
446  - (met_h_eta)*GradWz ) * mfy );
447  } else {
448  tau23(i,j,k) = myhalf * ( dv_dz
449  + ( (w(i, j, k) - w(i, j-1, k ))*dxInv[1]
450  - (met_h_eta)*GradWz ) * mfy );
451  }
452  tau32(i,j,k) = tau23(i,j,k);
453 
454  if (tau23i) tau23i(i,j,k) = myhalf * dv_dz;
455  });
456  }
457 
458  //***********************************************************************************
459  // Z-Dirichlet
460  //***********************************************************************************
461  if (zl_u_dir) {
462  Box planexz = tbxxz; planexz.setBig(2, planexz.smallEnd(2) );
463  tbxxz.growLo(2,-1);
464 
465  ParallelFor(planexz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
466  // Third order stencil with variable dz
467  Real dz0 = myhalf * ( z_nd(i,j,k+1) + z_nd(i,j+1,k+1)
468  - z_nd(i,j,k ) - z_nd(i,j+1,k ) );
469  Real dz1 = myhalf * ( z_nd(i,j,k+2) + z_nd(i,j+1,k+2)
470  - z_nd(i,j,k+1) - z_nd(i,j+1,k+1) );
471  Real idz0 = one / dz0;
472  Real f = (dz1 / dz0) + two;
473  Real f2 = f*f;
474  Real c3 = two / (f - f2);
475  Real c2 = -f2*c3;
476  Real c1 = -(one-f2)*c3;
477 
478  Real GradWz = myhalf * idz0 * ( w(i,j,k+1) + w(i-1,j,k+1)
479  - w(i,j,k ) - w(i-1,j,k ) );
480 
481  Real mfx = mf_ux(i,j,0);
482 
483  Real met_h_xi = Compute_h_xi_AtEdgeCenterJ (i,j,k,dxInv,z_nd);
484 
485  // Note: u(i,j,k) and u(i,j,k+1) are located at cell centers
486  // whereas u(i,j,k-1) is the value on the boundary
487  Real du_dz = (c1 * u(i,j,k-1) + c2 * u(i,j,k) + c3 * u(i,j,k+1))*idz0;
488  tau13(i,j,k) = myhalf * ( du_dz
489  + ( (w(i, j, k) - w(i-1, j, k))*dxInv[0]
490  - (met_h_xi)*GradWz ) * mfx );
491  tau31(i,j,k) = tau13(i,j,k);
492 
493  if (tau13i) tau13i(i,j,k) = myhalf * du_dz;
494  });
495  }
496  if (zh_u_dir) {
497  // NOTE: h_xi = 0
498  Box planexz = tbxxz; planexz.setSmall(2, planexz.bigEnd(2) );
499  tbxxz.growHi(2,-1);
500 
501  ParallelFor(planexz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
502  // Third order stencil with variable dz
503  Real dz0 = myhalf * ( z_nd(i,j,k ) + z_nd(i,j+1,k )
504  - z_nd(i,j,k-1) - z_nd(i,j+1,k-1) );
505  Real dz1 = myhalf * ( z_nd(i,j,k-1) + z_nd(i,j+1,k-1)
506  - z_nd(i,j,k-2) - z_nd(i,j+1,k-2) );
507  Real idz0 = one / dz0;
508  Real f = (dz1 / dz0) + two;
509  Real f2 = f*f;
510  Real c3 = two / (f - f2);
511  Real c2 = -f2*c3;
512  Real c1 = -(one-f2)*c3;
513 
514  Real mfx = mf_ux(i,j,0);
515 
516  Real du_dz = -(c1 * u(i,j,k) + c2 * u(i,j,k-1) + c3 * u(i,j,k-2))*idz0;
517  tau13(i,j,k) = myhalf * ( du_dz
518  + (w(i, j, k) - w(i-1, j, k))*dxInv[0]*mfx );
519  tau31(i,j,k) = tau13(i,j,k);
520 
521  if (tau13i) tau13i(i,j,k) = myhalf * du_dz;
522  });
523  }
524 
525  if (zl_v_dir) {
526  Box planeyz = tbxyz; planeyz.setBig(2, planeyz.smallEnd(2) );
527  tbxyz.growLo(2,-1);
528 
529  ParallelFor(planeyz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
530  // Third order stencil with variable dz
531  Real dz0 = myhalf * ( z_nd(i,j,k+1) + z_nd(i+1,j,k+1)
532  - z_nd(i,j,k ) - z_nd(i+1,j,k ) );
533  Real dz1 = myhalf * ( z_nd(i,j,k+2) + z_nd(i+1,j,k+2)
534  - z_nd(i,j,k+1) - z_nd(i+1,j,k+1) );
535  Real idz0 = one / dz0;
536  Real f = (dz1 / dz0) + two;
537  Real f2 = f*f;
538  Real c3 = two / (f - f2);
539  Real c2 = -f2*c3;
540  Real c1 = -(one-f2)*c3;
541 
542  Real GradWz = myhalf * idz0 * ( w(i ,j ,k+1) + w(i ,j-1,k+1)
543  - w(i ,j ,k ) - w(i ,j-1,k ) );
544 
545  Real mfy = mf_vy(i,j,0);
546 
547  Real met_h_eta;
548  met_h_eta = Compute_h_eta_AtEdgeCenterI (i,j,k,dxInv,z_nd);
549 
550  // Note: v(i,j,k) and v(i,j,k+1) are located at cell centers
551  // whereas v(i,j,k-1) is the value on the boundary
552  Real dv_dz = (c1 * v(i,j,k-1) + c2 * v(i,j,k ) + c3 * v(i,j,k+1))*idz0;
553  tau23(i,j,k) = myhalf * ( dv_dz
554  + ( (w(i, j, k) - w(i, j-1, k))*dxInv[1]
555  - (met_h_eta)*GradWz ) * mfy );
556  tau32(i,j,k) = tau23(i,j,k);
557 
558  if (tau23i) tau23i(i,j,k) = myhalf * dv_dz;
559  });
560  }
561  if (zh_v_dir && (tbxyz.bigEnd(2) == domain_yz.bigEnd(2))) {
562  // NOTE: h_eta = 0
563  Box planeyz = tbxyz; planeyz.setSmall(2, planeyz.bigEnd(2) );
564  tbxyz.growHi(2,-1);
565 
566  ParallelFor(planeyz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
567  // Third order stencil with variable dz
568  Real dz0 = myhalf * ( z_nd(i,j,k ) + z_nd(i+1,j,k )
569  - z_nd(i,j,k-1) - z_nd(i+1,j,k-1) );
570  Real dz1 = myhalf * ( z_nd(i,j,k-1) + z_nd(i+1,j,k-1)
571  - z_nd(i,j,k-2) - z_nd(i+1,j,k-2) );
572  Real idz0 = one / dz0;
573  Real f = (dz1 / dz0) + two;
574  Real f2 = f*f;
575  Real c3 = two / (f - f2);
576  Real c2 = -f2*c3;
577  Real c1 = -(one-f2)*c3;
578 
579  Real mfy = mf_vy(i,j,0);
580 
581  Real dv_dz = -(c1 * v(i,j,k ) + c2 * v(i,j,k-1) + c3 * v(i,j,k-2))*idz0;
582  tau23(i,j,k) = myhalf * ( dv_dz
583  + (w(i, j, k) - w(i, j-1, k))*dxInv[1]*mfy );
584  tau32(i,j,k) = tau23(i,j,k);
585 
586  if (tau23i) tau23i(i,j,k) = myhalf * dv_dz;
587  });
588  }
589 
590  // HO derivatives w/ Dirichlet BC (\partial <var> / \partial z from terrain transform)
591  if (zl_u_dir && zl_v_dir) {
592  Box planecc = bxcc; planecc.setBig(2, planecc.smallEnd(2) );
593  bxcc.growLo(2,-1);
594 
595  ParallelFor(planecc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
596  // Third order stencil with variable dz
597  Real dz0 = fourth * ( z_nd(i,j,k+1) + z_nd(i,j+1,k+1) + z_nd(i+1,j,k+1) + z_nd(i+1,j+1,k+1)
598  - z_nd(i,j,k ) - z_nd(i,j+1,k ) - z_nd(i+1,j,k ) - z_nd(i+1,j+1,k ) );
599  Real dz1 = fourth * ( z_nd(i,j,k+2) + z_nd(i,j+1,k+2) + z_nd(i+1,j,k+2) + z_nd(i+1,j+1,k+2)
600  - z_nd(i,j,k+1) - z_nd(i,j+1,k+1) - z_nd(i+1,j,k+1) - z_nd(i+1,j+1,k+1) );
601  Real idz0 = one / dz0;
602  Real f = (dz1 / dz0) + two;
603  Real f2 = f*f;
604  Real c3 = two / (f - f2);
605  Real c2 = -f2*c3;
606  Real c1 = -(one-f2)*c3;
607 
608  Real GradUz = myhalf * idz0 * ( (c1 * u(i ,j,k-1) + c2 * u(i ,j,k) + c3 * u(i ,j,k+1))
609  + (c1 * u(i-1,j,k-1) + c2 * u(i-1,j,k) + c3 * u(i-1,j,k+1)) );
610  Real GradVz = myhalf * idz0 * ( (c1 * v(i,j ,k-1) + c2 * v(i,j ,k) + c3 * v(i,j ,k+1))
611  + (c1 * v(i,j-1,k-1) + c2 * v(i,j-1,k) + c3 * v(i,j-1,k+1)) );
612 
613  Real mfx = mf_mx(i,j,0);
614  Real mfy = mf_my(i,j,0);
615 
616  Real met_h_xi,met_h_eta;
617  met_h_xi = Compute_h_xi_AtCellCenter (i,j,k,dxInv,z_nd);
618  met_h_eta = Compute_h_eta_AtCellCenter (i,j,k,dxInv,z_nd);
619 
620  tau11(i,j,k) = ( (u(i+1, j, k) - u(i, j, k) )*dxInv[0]
621  - (met_h_xi)*GradUz ) * mfx;
622  tau22(i,j,k) = ( (v(i, j+1, k) - v(i, j, k) )*dxInv[1]
623  - (met_h_eta)*GradVz ) * mfy;
624  tau33(i,j,k) = ( w(i, j, k+1) - w(i, j, k) )*idz0;
625 
626  });
627 
628  Box planexy = tbxxy; planexy.setBig(2, planexy.smallEnd(2) );
629  tbxxy.growLo(2,-1);
630 
631  ParallelFor(planexy,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
632  // Third order stencil with variable dz
633  Real dz0 = ( z_nd(i,j,k+1) - z_nd(i,j,k ) );
634  Real dz1 = ( z_nd(i,j,k+2) - z_nd(i,j,k+1) );
635  Real idz0 = one / dz0;
636  Real f = (dz1 / dz0) + two;
637  Real f2 = f*f;
638  Real c3 = two / (f - f2);
639  Real c2 = -f2*c3;
640  Real c1 = -(one-f2)*c3;
641 
642  Real GradUz = myhalf * idz0 * ( (c1 * u(i,j ,k-1) + c2 * u(i,j ,k) + c3 * u(i,j ,k+1))
643  + (c1 * u(i,j-1,k-1) + c2 * u(i,j-1,k) + c3 * u(i,j-1,k+1)) );
644  Real GradVz = myhalf * idz0 * ( (c1 * v(i ,j,k-1) + c2 * v(i ,j,k) + c3 * v(i ,j,k+1))
645  + (c1 * v(i-1,j,k-1) + c2 * v(i-1,j,k) + c3 * v(i-1,j,k+1)) );
646 
647  Real mfy = myhalf * (mf_uy(i,j,0) + mf_uy(i ,j-1,0));
648  Real mfx = myhalf * (mf_vx(i,j,0) + mf_vx(i-1,j ,0));
649 
650  Real met_h_xi,met_h_eta;
651  met_h_xi = Compute_h_xi_AtEdgeCenterK (i,j,k,dxInv,z_nd);
652  met_h_eta = Compute_h_eta_AtEdgeCenterK (i,j,k,dxInv,z_nd);
653 
654  tau12(i,j,k) = myhalf * ( (u(i, j, k) - u(i , j-1, k))*dxInv[1]*mfy
655  + (v(i, j, k) - v(i-1, j , k))*dxInv[0]*mfx
656  - (met_h_eta)*GradUz*mfy
657  - (met_h_xi )*GradVz*mfx );
658  tau21(i,j,k) = tau12(i,j,k);
659  });
660  }
661 
662  //***********************************************************************************
663  // Z-lo w/out Z-Dirichlet (GradWz extrapolation)
664  //***********************************************************************************
665  if (!zl_u_dir && (tbxxz.smallEnd(2) == domain_xz.smallEnd(2)) ) {
666  Box planexz = tbxxz; planexz.setBig(2, planexz.smallEnd(2));
667  tbxxz.growLo(2,-1);
668 
669  ParallelFor(planexz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
670  Real mfx = mf_ux(i,j,0);
671 
672  Real met_h_xi,met_h_zeta;
673  met_h_xi = Compute_h_xi_AtEdgeCenterJ (i,j,k,dxInv,z_nd);
674  met_h_zeta = Compute_h_zeta_AtEdgeCenterJ(i,j,k,dxInv,z_nd);
675 
676  Real GradWz = myhalf * dxInv[2] * ( w(i ,j ,k+1) + w(i-1,j ,k+1)
677  - w(i ,j ,k ) - w(i-1,j ,k ) );
678  GradWz /= met_h_zeta;
679 
680  Real du_dz = (u(i, j, k) - u(i , j, k-1))*dxInv[2]/met_h_zeta;
681  tau13(i,j,k) = myhalf * ( du_dz
682  + ( (w(i, j, k) - w(i-1, j, k ))*dxInv[0]
683  - (met_h_xi)*GradWz ) * mfx);
684  tau31(i,j,k) = tau13(i,j,k);
685 
686  if (tau13i) tau13i(i,j,k) = myhalf * du_dz;
687  });
688  }
689  if (!zl_v_dir && (tbxyz.smallEnd(2) == domain_yz.smallEnd(2))) {
690  Box planeyz = tbxyz; planeyz.setBig(2, planeyz.smallEnd(2) );
691  tbxyz.growLo(2,-1);
692 
693  ParallelFor(planeyz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
694  Real mfy = mf_vy(i,j,0);
695 
696  Real met_h_eta,met_h_zeta;
697  met_h_eta = Compute_h_eta_AtEdgeCenterI (i,j,k,dxInv,z_nd);
698  met_h_zeta = Compute_h_zeta_AtEdgeCenterI(i,j,k,dxInv,z_nd);
699 
700  Real GradWz = myhalf * dxInv[2] * ( w(i ,j ,k+1) + w(i ,j-1,k+1)
701  - w(i ,j ,k ) - w(i ,j-1,k ) );
702  GradWz /= met_h_zeta;
703 
704  Real dv_dz = (v(i, j, k) - v(i, j , k-1))*dxInv[2]/met_h_zeta;
705  tau23(i,j,k) = myhalf * ( dv_dz
706  + ( (w(i, j, k) - w(i, j-1, k ))*dxInv[1]
707  - (met_h_eta)*GradWz ) * mfy );
708  tau32(i,j,k) = tau23(i,j,k);
709 
710  if (tau23i) tau23i(i,j,k) = myhalf * dv_dz;
711  });
712  }
713 
714  //***********************************************************************************
715  // Z-hi w/out Z-Dirichlet (h_xi = h_eta = 0)
716  //***********************************************************************************
717  if (!zh_u_dir && (tbxxz.bigEnd(2) == domain_xz.bigEnd(2))) {
718  Box planexz = tbxxz; planexz.setSmall(2, planexz.bigEnd(2) );
719  tbxxz.growHi(2,-1);
720 
721  ParallelFor(planexz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
722  Real mfx = mf_ux(i,j,0);
723 
724  Real met_h_zeta;
725  met_h_zeta = Compute_h_zeta_AtEdgeCenterJ(i,j,k,dxInv,z_nd);
726 
727  Real du_dz = (u(i, j, k) - u(i , j, k-1))*dxInv[2]/met_h_zeta;
728  tau13(i,j,k) = myhalf * ( du_dz
729  + (w(i, j, k) - w(i-1, j, k ))*dxInv[0]*mfx );
730  tau31(i,j,k) = tau13(i,j,k);
731 
732  if (tau13i) tau13i(i,j,k) = myhalf * du_dz;
733  });
734  }
735  if (!zh_v_dir && (tbxyz.bigEnd(2) == domain_yz.bigEnd(2))) {
736  Box planeyz = tbxyz; planeyz.setSmall(2, planeyz.bigEnd(2) );
737  tbxyz.growHi(2,-1);
738 
739  ParallelFor(planeyz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
740  Real mfy = mf_vy(i,j,0);
741 
742  Real met_h_zeta;
743  met_h_zeta = Compute_h_zeta_AtEdgeCenterI(i,j,k,dxInv,z_nd);
744 
745  Real dv_dz = (v(i, j, k) - v(i, j , k-1))*dxInv[2]/met_h_zeta;
746  tau23(i,j,k) = myhalf * ( dv_dz
747  + (w(i, j, k) - w(i, j-1, k ))*dxInv[1]*mfy );
748  tau32(i,j,k) = tau23(i,j,k);
749 
750  if (tau23i) tau23i(i,j,k) = myhalf * dv_dz;
751  });
752  }
753 
754  //***********************************************************************************
755  // Fill the interior cells
756  //***********************************************************************************
757  // Cell centered strains
758  ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
759  Real dz0 = fourth * ( z_nd(i,j,k+1) + z_nd(i,j+1,k+1) + z_nd(i+1,j,k+1) + z_nd(i+1,j+1,k+1)
760  - z_nd(i,j,k-1) - z_nd(i,j+1,k-1) - z_nd(i+1,j,k-1) - z_nd(i+1,j+1,k-1) );
761  Real idz0 = one / dz0;
762 
763  Real GradUz = (k == 0) ?
764  idz0 * ( u(i ,j ,k+1) + u(i-1,j ,k+1)
765  - u(i ,j ,k ) - u(i-1,j ,k ) ) :
766  myhalf * idz0 * ( u(i ,j ,k+1) + u(i-1,j ,k+1)
767  - u(i ,j ,k-1) - u(i-1,j ,k-1) );
768  Real GradVz = (k == 0) ?
769  idz0 * ( v(i ,j ,k+1) + v(i ,j-1,k+1)
770  - v(i ,j ,k ) - v(i ,j-1,k ) ) :
771  myhalf * idz0 * ( v(i ,j ,k+1) + v(i ,j-1,k+1)
772  - v(i ,j ,k-1) - v(i ,j-1,k-1) );
773 
774  Real mfx = mf_mx(i,j,0);
775  Real mfy = mf_my(i,j,0);
776 
777  Real met_h_xi,met_h_eta,met_h_zeta;
778  met_h_xi = Compute_h_xi_AtCellCenter (i,j,k,dxInv,z_nd);
779  met_h_eta = Compute_h_eta_AtCellCenter (i,j,k,dxInv,z_nd);
780  met_h_zeta = detJ(i,j,k);
781 
782  tau11(i,j,k) = ( (u(i+1, j, k) - u(i, j, k))*dxInv[0] - met_h_xi*GradUz ) * mfx;
783  tau22(i,j,k) = ( (v(i, j+1, k) - v(i, j, k))*dxInv[1] - met_h_eta*GradVz ) * mfy;
784  tau33(i,j,k) = ( w(i, j, k+1) - w(i, j, k) )*dxInv[2]/met_h_zeta;
785  });
786 
787  // Off-diagonal strains
788  ParallelFor(tbxxy,tbxxz,tbxyz,
789  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
790  Real dz0 = ( z_nd(i,j,k+1) - z_nd(i,j,k-1) );
791  Real idz0 = one / dz0;
792 
793  Real GradUz = (k == 0) ?
794  idz0 * ( u(i ,j ,k+1) + u(i ,j-1,k+1)
795  - u(i ,j ,k ) - u(i ,j-1,k ) ) :
796  myhalf * idz0 * ( u(i ,j ,k+1) + u(i ,j-1,k+1)
797  - u(i ,j ,k-1) - u(i ,j-1,k-1) );
798  Real GradVz = (k == 0) ?
799  idz0 * ( v(i ,j ,k+1) + v(i-1,j ,k+1)
800  - v(i ,j ,k ) - v(i-1,j ,k ) ) :
801  myhalf * idz0 * ( v(i ,j ,k+1) + v(i-1,j ,k+1)
802  - v(i ,j ,k-1) - v(i-1,j ,k-1) );
803 
804  Real mfy = myhalf * (mf_uy(i,j,0) + mf_uy(i ,j-1,0));
805  Real mfx = myhalf * (mf_vx(i,j,0) + mf_vx(i-1,j ,0));
806 
807  Real met_h_xi,met_h_eta;
808  met_h_xi = Compute_h_xi_AtEdgeCenterK (i,j,k,dxInv,z_nd);
809  met_h_eta = Compute_h_eta_AtEdgeCenterK (i,j,k,dxInv,z_nd);
810 
811  tau12(i,j,k) = myhalf * ( (u(i, j, k) - u(i , j-1, k))*dxInv[1]*mfy
812  + (v(i, j, k) - v(i-1, j , k))*dxInv[0]*mfx
813  - (met_h_eta)*GradUz*mfy
814  - (met_h_xi )*GradVz*mfx );
815  tau21(i,j,k) = tau12(i,j,k);
816  },
817  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
818  Real dz0 = myhalf * ( z_nd(i,j,k+1) + z_nd(i,j+1,k+1)
819  - z_nd(i,j,k-1) - z_nd(i,j+1,k-1) );
820  Real idz0 = one / dz0;
821 
822  Real GradWz = myhalf * idz0 * ( w(i ,j ,k+1) + w(i-1,j ,k+1)
823  - w(i ,j ,k-1) - w(i-1,j ,k-1) );
824 
825  Real mfx = mf_ux(i,j,0);
826 
827  Real met_h_xi,met_h_zeta;
828  met_h_xi = Compute_h_xi_AtEdgeCenterJ (i,j,k,dxInv,z_nd);
829  met_h_zeta = Compute_h_zeta_AtEdgeCenterJ(i,j,k,dxInv,z_nd);
830 
831  Real du_dz = (u(i, j, k) - u(i , j, k-1))*dxInv[2]/met_h_zeta;
832  tau13(i,j,k) = myhalf * ( du_dz
833  + ( (w(i, j, k) - w(i-1, j, k ))*dxInv[0]
834  - (met_h_xi)*GradWz ) * mfx );
835  tau31(i,j,k) = tau13(i,j,k);
836 
837  if (tau13i) tau13i(i,j,k) = myhalf * du_dz;
838  },
839  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
840  Real dz0 = myhalf * ( z_nd(i,j,k+1) + z_nd(i+1,j,k+1)
841  - z_nd(i,j,k-1) - z_nd(i+1,j,k-1) );
842  Real idz0 = one / dz0;
843 
844  Real GradWz = myhalf * idz0 * ( w(i ,j ,k+1) + w(i ,j-1,k+1)
845  - w(i ,j ,k-1) - w(i ,j-1,k-1) );
846 
847  Real mfy = mf_vy(i,j,0);
848 
849  Real met_h_eta,met_h_zeta;
850  met_h_eta = Compute_h_eta_AtEdgeCenterI (i,j,k,dxInv,z_nd);
851  met_h_zeta = Compute_h_zeta_AtEdgeCenterI(i,j,k,dxInv,z_nd);
852 
853  Real dv_dz = (v(i, j, k) - v(i, j , k-1))*dxInv[2]/met_h_zeta;
854  tau23(i,j,k) = myhalf * ( dv_dz
855  + ( (w(i, j, k) - w(i, j-1, k ))*dxInv[1]
856  - (met_h_eta)*GradWz ) * mfy );
857  tau32(i,j,k) = tau23(i,j,k);
858 
859  if (tau23i) tau23i(i,j,k) = myhalf * dv_dz;
860  });
861 }
three
constexpr amrex::Real three
Definition: ERF_Constants.H:11
two
constexpr amrex::Real two
Definition: ERF_Constants.H:10
one
constexpr amrex::Real one
Definition: ERF_Constants.H:9
fourth
constexpr amrex::Real fourth
Definition: ERF_Constants.H:14
zero
constexpr amrex::Real zero
Definition: ERF_Constants.H:8
myhalf
constexpr amrex::Real myhalf
Definition: ERF_Constants.H:13
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
tau32
@ tau32
Definition: ERF_DataStruct.H:38
tau31
@ tau31
Definition: ERF_DataStruct.H:38
tau21
@ tau21
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
Compute_h_xi_AtEdgeCenterK
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_xi_AtEdgeCenterK(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:243
Compute_h_xi_AtEdgeCenterJ
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_xi_AtEdgeCenterJ(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:289
Compute_h_zeta_AtEdgeCenterJ
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_zeta_AtEdgeCenterJ(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:275
Compute_h_zeta_AtEdgeCenterI
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_zeta_AtEdgeCenterI(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:319
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:85
Compute_h_eta_AtEdgeCenterK
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_eta_AtEdgeCenterK(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:258
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:70
Compute_h_eta_AtEdgeCenterI
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_eta_AtEdgeCenterI(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:346
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
module_model_constants::c2
real(c_double), parameter c2
Definition: ERF_module_model_constants.F90:35
module_mp_morr_two_moment::c1
real(c_double), private c1
Definition: ERF_module_mp_morr_two_moment.F90:212

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: