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