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 bx, 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 > &SmnSmn_a, const Real)
 

Function Documentation

◆ ComputeStrain_T()

void ComputeStrain_T ( Box  bx,
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 > &  SmnSmn_a,
const  Real 
)

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]implicit_fac– factor of implicitness for vertical differences only
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  } else {
163  tau12(i,j,k) = 0.5 * ( (u(i, j, k) - u(i , j-1, k))*dxInv[1]*mfy
164  + (v(i, j, k) - v(i-1, j , k))*dxInv[0]*mfx
165  - (met_h_eta)*GradUz*mfy
166  - (met_h_xi )*GradVz*mfx );
167  }
168  tau21(i,j,k) = tau12(i,j,k);
169  });
170  }
171  if (xh_v_dir) {
172  Box planexy = tbxxy; planexy.setSmall(0, planexy.bigEnd(0) );
173  tbxxy.growHi(0,-1);
174  bool need_to_test = (bc_ptr[BCVars::yvel_bc].hi(0) == ERFBCType::ext_dir_upwind) ? true : false;
175 
176  ParallelFor(planexy,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
177  Real inv_dist = (k == 0) ? Real(2.0) / (z_nd(i,j,k+2) - z_nd(i,j,k )) :
178  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));
179 
180  Real GradUz = (k == 0) ?
181  inv_dist * ( u(i ,j ,k+1) + u(i ,j-1,k+1)
182  - u(i ,j ,k ) - u(i ,j-1,k ) ) :
183  inv_dist * ( u(i ,j ,k+1) + u(i ,j-1,k+1)
184  - u(i ,j ,k-1) - u(i ,j-1,k-1) );
185  Real GradVz = (k == 0) ?
186  inv_dist * ( v(i ,j ,k+1) + v(i-1,j ,k+1)
187  - v(i ,j ,k ) - v(i-1,j ,k ) ) :
188  inv_dist * ( v(i ,j ,k+1) + v(i-1,j ,k+1)
189  - v(i ,j ,k-1) - v(i-1,j ,k-1) );
190 
191  Real mfy = 0.5 * (mf_uy(i,j,0) + mf_uy(i ,j-1,0));
192  Real mfx = 0.5 * (mf_vx(i,j,0) + mf_vx(i-1,j ,0));
193 
194  Real met_h_xi = Compute_h_xi_AtEdgeCenterK (i,j,k,dxInv,z_nd);
195  Real met_h_eta = Compute_h_eta_AtEdgeCenterK (i,j,k,dxInv,z_nd);
196 
197  if (!need_to_test || u(dom_hi.x+1,j,k) <= 0.) {
198  tau12(i,j,k) = 0.5 * ( (u(i, j, k) - u(i, j-1, k) )*dxInv[1]*mfy
199  - (-(8./3.) * v(i,j,k) + 3. * v(i-1,j,k) - (1./3.) * v(i-2,j,k))*dxInv[0]*mfx
200  - (met_h_eta)*GradUz*mfy
201  - (met_h_xi )*GradVz*mfx );
202  } else {
203  tau12(i,j,k) = 0.5 * ( (u(i, j, k) - u(i , j-1, k))*dxInv[1]*mfy
204  + (v(i, j, k) - v(i-1, j , k))*dxInv[0]*mfx
205  - (met_h_eta)*GradUz*mfy
206  - (met_h_xi )*GradVz*mfx );
207  }
208  tau21(i,j,k) = tau12(i,j,k);
209  });
210  }
211 
212  if (xl_w_dir) {
213  Box planexz = tbxxz; planexz.setBig(0, planexz.smallEnd(0) );
214  planexz.setSmall(2, planexz.smallEnd(2)+1 ); planexz.setBig(2, planexz.bigEnd(2)-1 );
215  tbxxz.growLo(0,-1);
216  bool need_to_test = (bc_ptr[BCVars::zvel_bc].lo(0) == ERFBCType::ext_dir_upwind) ? true : false;
217 
218  ParallelFor(planexz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
219  Real dz0 = 0.5 * ( z_nd(i,j,k+1) + z_nd(i,j+1,k+1)
220  - z_nd(i,j,k-1) - z_nd(i,j+1,k-1) );
221  Real idz0 = 1.0 / dz0;
222 
223  Real GradWz = 0.5 * idz0 * ( w(i ,j ,k+1) + w(i-1,j ,k+1)
224  - w(i ,j ,k-1) - w(i-1,j ,k-1) );
225  Real mfx = mf_ux(i,j,0);
226 
227  Real met_h_xi = Compute_h_xi_AtEdgeCenterJ (i,j,k,dxInv,z_nd);
228  Real met_h_zeta = Compute_h_zeta_AtEdgeCenterJ(i,j,k,dxInv,z_nd);
229 
230  if (!need_to_test || u(dom_lo.x,j,k) <= 0.) {
231  tau13(i,j,k) = 0.5 * ( (u(i, j, k) - u(i, j, k-1))*dxInv[2]/met_h_zeta
232  + ( (-(8./3.) * w(i-1,j,k) + 3. * w(i,j,k) - (1./3.) * w(i+1,j,k))*dxInv[0]
233  - (met_h_xi)*GradWz ) * mfx );
234  } else {
235  tau13(i,j,k) = 0.5 * ( (u(i, j, k) - u(i , j, k-1))*dxInv[2]/met_h_zeta
236  + ( (w(i, j, k) - w(i-1, j, k ))*dxInv[0]
237  - (met_h_xi)*GradWz ) * mfx );
238  }
239  tau31(i,j,k) = tau13(i,j,k);
240  });
241  }
242 
243  if (xh_w_dir) {
244  Box planexz = tbxxz; planexz.setSmall(0, planexz.bigEnd(0) );
245  planexz.setSmall(2, planexz.smallEnd(2)+1 ); planexz.setBig(2, planexz.bigEnd(2)-1 );
246  tbxxz.growHi(0,-1);
247  bool need_to_test = (bc_ptr[BCVars::zvel_bc].hi(0) == ERFBCType::ext_dir_upwind) ? true : false;
248 
249  ParallelFor(planexz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
250  Real dz0 = 0.5 * ( z_nd(i,j,k+1) + z_nd(i,j+1,k+1)
251  - z_nd(i,j,k-1) - z_nd(i,j+1,k-1) );
252  Real idz0 = 1.0 / dz0;
253 
254  Real GradWz = 0.5 * idz0 * ( w(i ,j ,k+1) + w(i-1,j ,k+1)
255  - w(i ,j ,k-1) - w(i-1,j ,k-1) );
256 
257  Real mfx = mf_ux(i,j,0);
258 
259  Real met_h_xi = Compute_h_xi_AtEdgeCenterJ (i,j,k,dxInv,z_nd);
260  Real met_h_zeta = Compute_h_zeta_AtEdgeCenterJ(i,j,k,dxInv,z_nd);
261 
262  if (!need_to_test || u(dom_hi.x+1,j,k) <= 0.) {
263  tau13(i,j,k) = 0.5 * ( (u(i, j, k) - u(i, j, k-1))*dxInv[2]/met_h_zeta
264  - ( (-(8./3.) * w(i,j,k) + 3. * w(i-1,j,k) - (1./3.) * w(i-2,j,k))*dxInv[0]
265  - (met_h_xi)*GradWz ) * mfx );
266  } else {
267  tau13(i,j,k) = 0.5 * ( (u(i, j, k) - u(i , j, k-1))*dxInv[2]/met_h_zeta
268  + ( (w(i, j, k) - w(i-1, j, k ))*dxInv[0]
269  - (met_h_xi)*GradWz ) * mfx );
270  }
271  tau31(i,j,k) = tau13(i,j,k);
272  });
273  }
274 
275  //***********************************************************************************
276  // Y-Dirichlet
277  //***********************************************************************************
278  if (yl_u_dir) {
279  Box planexy = tbxxy; planexy.setBig(1, planexy.smallEnd(1) );
280  tbxxy.growLo(1,-1);
281  bool need_to_test = (bc_ptr[BCVars::xvel_bc].lo(1) == ERFBCType::ext_dir_upwind) ? true : false;
282 
283  ParallelFor(planexy,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
284  Real dz0 = ( z_nd(i,j,k+1) - z_nd(i,j,k-1) );
285  Real idz0 = 1.0 / dz0;
286 
287  Real GradUz = (k == 0) ?
288  idz0 * ( u(i ,j ,k+1) + u(i ,j-1,k+1)
289  - u(i ,j ,k ) - u(i ,j-1,k ) ) :
290  0.5 * idz0 * ( u(i ,j ,k+1) + u(i ,j-1,k+1)
291  - u(i ,j ,k-1) - u(i ,j-1,k-1) );
292  Real GradVz = (k == 0) ?
293  idz0 * ( v(i ,j ,k+1) + v(i-1,j ,k+1)
294  - v(i ,j ,k ) - v(i-1,j ,k ) ) :
295  0.5 * idz0 * ( v(i ,j ,k+1) + v(i-1,j ,k+1)
296  - v(i ,j ,k-1) - v(i-1,j ,k-1) );
297 
298  Real mfy = 0.5 * (mf_uy(i,j,0) + mf_uy(i ,j-1,0));
299  Real mfx = 0.5 * (mf_vx(i,j,0) + mf_vx(i-1,j ,0));
300 
301  Real met_h_xi = Compute_h_xi_AtEdgeCenterK (i,j,k,dxInv,z_nd);
302  Real met_h_eta = Compute_h_eta_AtEdgeCenterK (i,j,k,dxInv,z_nd);
303 
304  if (!need_to_test || v(i,dom_lo.y,k) >= 0.) {
305  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
306  + (v(i, j, k) - v(i-1, j, k))*dxInv[0]*mfx
307  - (met_h_eta)*GradUz*mfy
308  - (met_h_xi )*GradVz*mfx );
309  } else {
310  tau12(i,j,k) = 0.5 * ( (u(i, j, k) - u(i , j-1, k))*dxInv[1]*mfy
311  + (v(i, j, k) - v(i-1, j , k))*dxInv[0]*mfx
312  - (met_h_eta)*GradUz*mfy
313  - (met_h_xi )*GradVz*mfx );
314  }
315  tau21(i,j,k) = tau12(i,j,k);
316  });
317  }
318  if (yh_u_dir) {
319  Box planexy = tbxxy; planexy.setSmall(1, planexy.bigEnd(1) );
320  tbxxy.growHi(1,-1);
321  bool need_to_test = (bc_ptr[BCVars::xvel_bc].hi(1) == ERFBCType::ext_dir_upwind) ? true : false;
322 
323  ParallelFor(planexy,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
324  Real dz0 = ( z_nd(i,j,k+1) - z_nd(i,j,k-1) );
325  Real idz0 = 1.0 / dz0;
326 
327  Real GradUz = (k == 0) ?
328  idz0 * ( u(i ,j ,k+1) + u(i ,j-1,k+1)
329  - u(i ,j ,k ) - u(i ,j-1,k ) ) :
330  0.5 * idz0 * ( u(i ,j ,k+1) + u(i ,j-1,k+1)
331  - u(i ,j ,k-1) - u(i ,j-1,k-1) );
332  Real GradVz = (k == 0) ?
333  idz0 * ( v(i ,j ,k+1) + v(i-1,j ,k+1)
334  - v(i ,j ,k ) - v(i-1,j ,k ) ) :
335  0.5 * idz0 * ( v(i ,j ,k+1) + v(i-1,j ,k+1)
336  - v(i ,j ,k-1) - v(i-1,j ,k-1) );
337 
338  Real mfy = 0.5 * (mf_uy(i,j,0) + mf_uy(i ,j-1,0));
339  Real mfx = 0.5 * (mf_vx(i,j,0) + mf_vx(i-1,j ,0));
340 
341  Real met_h_xi = Compute_h_xi_AtEdgeCenterK (i,j,k,dxInv,z_nd);
342  Real met_h_eta = Compute_h_eta_AtEdgeCenterK (i,j,k,dxInv,z_nd);
343 
344  if (!need_to_test || v(i,dom_hi.y+1,k) >= 0.) {
345  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 +
346  + (v(i, j, k) - v(i-1, j, k))*dxInv[0]*mfx
347  - (met_h_eta)*GradUz*mfy
348  - (met_h_xi )*GradVz*mfx );
349  } else {
350  tau12(i,j,k) = 0.5 * ( (u(i, j, k) - u(i , j-1, k))*dxInv[1]*mfy
351  + (v(i, j, k) - v(i-1, j , k))*dxInv[0]*mfx
352  - (met_h_eta)*GradUz*mfy
353  - (met_h_xi )*GradVz*mfx );
354  }
355  tau21(i,j,k) = tau12(i,j,k);
356  });
357  }
358 
359  if (yl_w_dir) {
360  Box planeyz = tbxyz; planeyz.setBig(1, planeyz.smallEnd(1) );
361  planeyz.setSmall(2, planeyz.smallEnd(2)+1 ); planeyz.setBig(2, planeyz.bigEnd(2)-1 );
362  tbxyz.growLo(1,-1);
363  bool need_to_test = (bc_ptr[BCVars::zvel_bc].lo(1) == ERFBCType::ext_dir_upwind) ? true : false;
364 
365  ParallelFor(planeyz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
366  Real dz0 = 0.5 * ( z_nd(i,j,k+1) + z_nd(i+1,j,k+1)
367  - z_nd(i,j,k-1) - z_nd(i+1,j,k-1) );
368  Real idz0 = 1.0 / dz0;
369 
370  Real GradWz = 0.5 * idz0 * ( w(i ,j ,k+1) + w(i ,j-1,k+1)
371  - w(i ,j ,k-1) - w(i ,j-1,k-1) );
372 
373  Real mfy = mf_vy(i,j,0);
374 
375  Real met_h_eta = Compute_h_eta_AtEdgeCenterI (i,j,k,dxInv,z_nd);
376  Real met_h_zeta = Compute_h_zeta_AtEdgeCenterI(i,j,k,dxInv,z_nd);
377 
378  if (!need_to_test || v(i,dom_lo.y,k) >= 0.) {
379  tau23(i,j,k) = 0.5 * ( (v(i, j, k) - v(i, j, k-1))*dxInv[2]/met_h_zeta
380  + ( (-(8./3.) * w(i,j-1,k) + 3. * w(i,j ,k) - (1./3.) * w(i,j+1,k))*dxInv[1]
381  - (met_h_eta)*GradWz ) * mfy );
382  } else {
383  tau23(i,j,k) = 0.5 * ( (v(i, j, k) - v(i, j , k-1))*dxInv[2]/met_h_zeta
384  + ( (w(i, j, k) - w(i, j-1, k ))*dxInv[1]
385  - (met_h_eta)*GradWz ) * mfy );
386  }
387  tau32(i,j,k) = tau23(i,j,k);
388  });
389  }
390  if (yh_w_dir) {
391  Box planeyz = tbxyz; planeyz.setSmall(1, planeyz.bigEnd(1) );
392  planeyz.setSmall(2, planeyz.smallEnd(2)+1 ); planeyz.setBig(2, planeyz.bigEnd(2)-1 );
393  tbxyz.growHi(1,-1);
394  bool need_to_test = (bc_ptr[BCVars::zvel_bc].hi(1) == ERFBCType::ext_dir_upwind) ? true : false;
395 
396  ParallelFor(planeyz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
397  Real dz0 = 0.5 * ( z_nd(i,j,k+1) + z_nd(i+1,j,k+1)
398  - z_nd(i,j,k-1) - z_nd(i+1,j,k-1) );
399  Real idz0 = 1.0 / dz0;
400 
401  Real GradWz = 0.5 * idz0 * ( w(i ,j ,k+1) + w(i ,j-1,k+1)
402  - w(i ,j ,k-1) - w(i ,j-1,k-1) );
403 
404  Real mfy = mf_vy(i,j,0);
405 
406  Real met_h_eta = Compute_h_eta_AtEdgeCenterI (i,j,k,dxInv,z_nd);
407  Real met_h_zeta = Compute_h_zeta_AtEdgeCenterI(i,j,k,dxInv,z_nd);
408 
409  if (!need_to_test || v(i,dom_hi.y+1,k) >= 0.) {
410  tau23(i,j,k) = 0.5 * ( (v(i, j, k) - v(i, j, k-1))*dxInv[2]/met_h_zeta
411  - ( (-(8./3.) * w(i,j ,k) + 3. * w(i,j-1,k) - (1./3.) * w(i,j-2,k))*dxInv[1]
412  - (met_h_eta)*GradWz ) * mfy );
413  } else {
414  tau23(i,j,k) = 0.5 * ( (v(i, j, k) - v(i, j , k-1))*dxInv[2]/met_h_zeta
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  }
421 
422  //***********************************************************************************
423  // Z-Dirichlet
424  //***********************************************************************************
425  if (zl_u_dir) {
426  Box planexz = tbxxz; planexz.setBig(2, planexz.smallEnd(2) );
427  tbxxz.growLo(2,-1);
428 
429  ParallelFor(planexz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
430  // Third order stencil with variable dz
431  Real dz0 = 0.5 * ( z_nd(i,j,k+1) + z_nd(i,j+1,k+1)
432  - z_nd(i,j,k ) - z_nd(i,j+1,k ) );
433  Real dz1 = 0.5 * ( z_nd(i,j,k+2) + z_nd(i,j+1,k+2)
434  - z_nd(i,j,k+1) - z_nd(i,j+1,k+1) );
435  Real idz0 = 1.0 / dz0;
436  Real f = (dz1 / dz0) + 2.0;
437  Real f2 = f*f;
438  Real c3 = 2.0 / (f - f2);
439  Real c2 = -f2*c3;
440  Real c1 = -(1.0-f2)*c3;
441 
442  Real GradWz = 0.5 * idz0 * ( w(i,j,k+1) + w(i-1,j,k+1)
443  - w(i,j,k ) - w(i-1,j,k ) );
444 
445  Real mfx = mf_ux(i,j,0);
446 
447  Real met_h_xi;
448  met_h_xi = Compute_h_xi_AtEdgeCenterJ (i,j,k,dxInv,z_nd);
449 
450  // Note: u(i,j,k) and u(i,j,k+1) are located at cell centers
451  // whereas u(i,j,k-1) is the value on the boundary
452  tau13(i,j,k) = 0.5 * ( (c1 * u(i,j,k-1) + c2 * u(i,j,k) + c3 * u(i,j,k+1))*idz0
453  + ( (w(i, j, k) - w(i-1, j, k))*dxInv[0]
454  - (met_h_xi)*GradWz ) * mfx );
455  tau31(i,j,k) = tau13(i,j,k);
456  });
457  }
458  if (zh_u_dir) {
459  // NOTE: h_xi = 0
460  Box planexz = tbxxz; planexz.setSmall(2, planexz.bigEnd(2) );
461  tbxxz.growHi(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 ) + z_nd(i,j+1,k )
466  - z_nd(i,j,k-1) - z_nd(i,j+1,k-1) );
467  Real dz1 = 0.5 * ( z_nd(i,j,k-1) + z_nd(i,j+1,k-1)
468  - z_nd(i,j,k-2) - z_nd(i,j+1,k-2) );
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 mfx = mf_ux(i,j,0);
477 
478  tau13(i,j,k) = 0.5 * ( -(c1 * u(i,j,k) + c2 * u(i,j,k-1) + c3 * u(i,j,k-2))*idz0
479  + (w(i, j, k) - w(i-1, j, k))*dxInv[0]*mfx );
480  tau31(i,j,k) = tau13(i,j,k);
481  });
482  }
483 
484  if (zl_v_dir) {
485  Box planeyz = tbxyz; planeyz.setBig(2, planeyz.smallEnd(2) );
486  tbxyz.growLo(2,-1);
487 
488  ParallelFor(planeyz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
489  // Third order stencil with variable dz
490  Real dz0 = 0.5 * ( z_nd(i,j,k+1) + z_nd(i+1,j,k+1)
491  - z_nd(i,j,k ) - z_nd(i+1,j,k ) );
492  Real dz1 = 0.5 * ( z_nd(i,j,k+2) + z_nd(i+1,j,k+2)
493  - z_nd(i,j,k+1) - z_nd(i+1,j,k+1) );
494  Real idz0 = 1.0 / dz0;
495  Real f = (dz1 / dz0) + 2.0;
496  Real f2 = f*f;
497  Real c3 = 2.0 / (f - f2);
498  Real c2 = -f2*c3;
499  Real c1 = -(1.0-f2)*c3;
500 
501  Real GradWz = 0.5 * idz0 * ( w(i ,j ,k+1) + w(i ,j-1,k+1)
502  - w(i ,j ,k ) - w(i ,j-1,k ) );
503 
504  Real mfy = mf_vy(i,j,0);
505 
506  Real met_h_eta;
507  met_h_eta = Compute_h_eta_AtEdgeCenterI (i,j,k,dxInv,z_nd);
508 
509  // Note: v(i,j,k) and v(i,j,k+1) are located at cell centers
510  // whereas v(i,j,k-1) is the value on the boundary
511  tau23(i,j,k) = 0.5 * ( (c1 * v(i,j,k-1) + c2 * v(i,j,k ) + c3 * v(i,j,k+1))*idz0
512  + ( (w(i, j, k) - w(i, j-1, k))*dxInv[1]
513  - (met_h_eta)*GradWz ) * mfy );
514  tau32(i,j,k) = tau23(i,j,k);
515  });
516  }
517  if (zh_v_dir && (tbxyz.bigEnd(2) == domain_yz.bigEnd(2))) {
518  // NOTE: h_eta = 0
519  Box planeyz = tbxyz; planeyz.setSmall(2, planeyz.bigEnd(2) );
520  tbxyz.growHi(2,-1);
521 
522  ParallelFor(planeyz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
523  // Third order stencil with variable dz
524  Real dz0 = 0.5 * ( z_nd(i,j,k ) + z_nd(i+1,j,k )
525  - z_nd(i,j,k-1) - z_nd(i+1,j,k-1) );
526  Real dz1 = 0.5 * ( z_nd(i,j,k-1) + z_nd(i+1,j,k-1)
527  - z_nd(i,j,k-2) - z_nd(i+1,j,k-2) );
528  Real idz0 = 1.0 / dz0;
529  Real f = (dz1 / dz0) + 2.0;
530  Real f2 = f*f;
531  Real c3 = 2.0 / (f - f2);
532  Real c2 = -f2*c3;
533  Real c1 = -(1.0-f2)*c3;
534 
535  Real mfy = mf_vy(i,j,0);
536 
537  tau23(i,j,k) = 0.5 * ( -(c1 * v(i,j,k ) + c2 * v(i,j,k-1) + c3 * v(i,j,k-2))*idz0
538  + (w(i, j, k) - w(i, j-1, k))*dxInv[1]*mfy );
539  tau32(i,j,k) = tau23(i,j,k);
540  });
541  }
542 
543  // HO derivatives w/ Dirichlet BC (\partival <var> / \partial z from terrain transform)
544  if (zl_u_dir && zl_v_dir) {
545  Box planecc = bxcc; planecc.setBig(2, planecc.smallEnd(2) );
546  bxcc.growLo(2,-1);
547 
548  ParallelFor(planecc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
549  // Third order stencil with variable dz
550  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)
551  - z_nd(i,j,k ) - z_nd(i,j+1,k ) - z_nd(i+1,j,k ) - z_nd(i+1,j+1,k ) );
552  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)
553  - 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) );
554  Real idz0 = 1.0 / dz0;
555  Real f = (dz1 / dz0) + 2.0;
556  Real f2 = f*f;
557  Real c3 = 2.0 / (f - f2);
558  Real c2 = -f2*c3;
559  Real c1 = -(1.0-f2)*c3;
560 
561  Real GradUz = 0.5 * idz0 * ( (c1 * u(i ,j,k-1) + c2 * u(i ,j,k) + c3 * u(i ,j,k+1))
562  + (c1 * u(i-1,j,k-1) + c2 * u(i-1,j,k) + c3 * u(i-1,j,k+1)) );
563  Real GradVz = 0.5 * idz0 * ( (c1 * v(i,j ,k-1) + c2 * v(i,j ,k) + c3 * v(i,j ,k+1))
564  + (c1 * v(i,j-1,k-1) + c2 * v(i,j-1,k) + c3 * v(i,j-1,k+1)) );
565 
566  Real mfx = mf_mx(i,j,0);
567  Real mfy = mf_my(i,j,0);
568 
569  Real met_h_xi,met_h_eta;
570  met_h_xi = Compute_h_xi_AtCellCenter (i,j,k,dxInv,z_nd);
571  met_h_eta = Compute_h_eta_AtCellCenter (i,j,k,dxInv,z_nd);
572 
573  tau11(i,j,k) = ( (u(i+1, j, k) - u(i, j, k) )*dxInv[0]
574  - (met_h_xi)*GradUz ) * mfx;
575  tau22(i,j,k) = ( (v(i, j+1, k) - v(i, j, k) )*dxInv[1]
576  - (met_h_eta)*GradVz ) * mfy;
577  tau33(i,j,k) = ( w(i, j, k+1) - w(i, j, k) )*idz0;
578 
579  });
580 
581  Box planexy = tbxxy; planexy.setBig(2, planexy.smallEnd(2) );
582  tbxxy.growLo(2,-1);
583 
584  ParallelFor(planexy,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
585  // Third order stencil with variable dz
586  Real dz0 = ( z_nd(i,j,k+1) - z_nd(i,j,k ) );
587  Real dz1 = ( z_nd(i,j,k+2) - z_nd(i,j,k+1) );
588  Real idz0 = 1.0 / dz0;
589  Real f = (dz1 / dz0) + 2.0;
590  Real f2 = f*f;
591  Real c3 = 2.0 / (f - f2);
592  Real c2 = -f2*c3;
593  Real c1 = -(1.0-f2)*c3;
594 
595  Real GradUz = 0.5 * idz0 * ( (c1 * u(i,j ,k-1) + c2 * u(i,j ,k) + c3 * u(i,j ,k+1))
596  + (c1 * u(i,j-1,k-1) + c2 * u(i,j-1,k) + c3 * u(i,j-1,k+1)) );
597  Real GradVz = 0.5 * idz0 * ( (c1 * v(i ,j,k-1) + c2 * v(i ,j,k) + c3 * v(i ,j,k+1))
598  + (c1 * v(i-1,j,k-1) + c2 * v(i-1,j,k) + c3 * v(i-1,j,k+1)) );
599 
600  Real mfy = 0.5 * (mf_uy(i,j,0) + mf_uy(i ,j-1,0));
601  Real mfx = 0.5 * (mf_vx(i,j,0) + mf_vx(i-1,j ,0));
602 
603  Real met_h_xi,met_h_eta;
604  met_h_xi = Compute_h_xi_AtEdgeCenterK (i,j,k,dxInv,z_nd);
605  met_h_eta = Compute_h_eta_AtEdgeCenterK (i,j,k,dxInv,z_nd);
606 
607  tau12(i,j,k) = 0.5 * ( (u(i, j, k) - u(i , j-1, k))*dxInv[1]*mfy
608  + (v(i, j, k) - v(i-1, j , k))*dxInv[0]*mfx
609  - (met_h_eta)*GradUz*mfy
610  - (met_h_xi )*GradVz*mfx );
611  tau21(i,j,k) = tau12(i,j,k);
612  });
613  }
614 
615  //***********************************************************************************
616  // Z-lo w/out Z-Dirichlet (GradWz extrapolation)
617  //***********************************************************************************
618  if (!zl_u_dir && (tbxxz.smallEnd(2) == domain_xz.smallEnd(2)) ) {
619  Box planexz = tbxxz; planexz.setBig(2, planexz.smallEnd(2));
620  tbxxz.growLo(2,-1);
621 
622  ParallelFor(planexz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
623  Real GradWz = 0.5 * dxInv[2] * ( w(i ,j ,k+1) + w(i-1,j ,k+1)
624  - w(i ,j ,k ) - w(i-1,j ,k ) );
625 
626  Real mfx = mf_ux(i,j,0);
627 
628  Real met_h_xi,met_h_zeta;
629  met_h_xi = Compute_h_xi_AtEdgeCenterJ (i,j,k,dxInv,z_nd);
630  met_h_zeta = Compute_h_zeta_AtEdgeCenterJ(i,j,k,dxInv,z_nd);
631 
632  tau13(i,j,k) = 0.5 * ( (u(i, j, k) - u(i , j, k-1))*dxInv[2]/met_h_zeta
633  + ( (w(i, j, k) - w(i-1, j, k ))*dxInv[0]
634  - (met_h_xi/met_h_zeta)*GradWz ) * mfx);
635  tau31(i,j,k) = tau13(i,j,k);
636  });
637  }
638  if (!zl_v_dir && (tbxyz.smallEnd(2) == domain_yz.smallEnd(2))) {
639  Box planeyz = tbxyz; planeyz.setBig(2, planeyz.smallEnd(2) );
640  tbxyz.growLo(2,-1);
641 
642  ParallelFor(planeyz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
643  Real GradWz = 0.5 * dxInv[2] * ( w(i ,j ,k+1) + w(i ,j-1,k+1)
644  - w(i ,j ,k ) - w(i ,j-1,k ) );
645 
646  Real mfy = mf_vy(i,j,0);
647 
648  Real met_h_eta,met_h_zeta;
649  met_h_eta = Compute_h_eta_AtEdgeCenterI (i,j,k,dxInv,z_nd);
650  met_h_zeta = Compute_h_zeta_AtEdgeCenterI(i,j,k,dxInv,z_nd);
651 
652  tau23(i,j,k) = 0.5 * ( (v(i, j, k) - v(i, j , k-1))*dxInv[2]/met_h_zeta
653  + ( (w(i, j, k) - w(i, j-1, k ))*dxInv[1]
654  - (met_h_eta/met_h_zeta)*GradWz ) * mfy );
655  tau32(i,j,k) = tau23(i,j,k);
656  });
657  }
658 
659  //***********************************************************************************
660  // Z-hi w/out Z-Dirichlet (h_xi = h_eta = 0)
661  //***********************************************************************************
662  if (!zh_u_dir && (tbxxz.bigEnd(2) == domain_xz.bigEnd(2))) {
663  Box planexz = tbxxz; planexz.setSmall(2, planexz.bigEnd(2) );
664  tbxxz.growHi(2,-1);
665 
666  ParallelFor(planexz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
667  Real mfx = mf_ux(i,j,0);
668 
669  Real met_h_zeta;
670  met_h_zeta = Compute_h_zeta_AtEdgeCenterJ(i,j,k,dxInv,z_nd);
671 
672  tau13(i,j,k) = 0.5 * ( (u(i, j, k) - u(i , j, k-1))*dxInv[2]/met_h_zeta
673  + (w(i, j, k) - w(i-1, j, k ))*dxInv[0]*mfx );
674  tau31(i,j,k) = tau13(i,j,k);
675  });
676  }
677  if (!zh_v_dir && (tbxyz.bigEnd(2) == domain_yz.bigEnd(2))) {
678  Box planeyz = tbxyz; planeyz.setSmall(2, planeyz.bigEnd(2) );
679  tbxyz.growHi(2,-1);
680 
681  ParallelFor(planeyz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
682  Real mfy = mf_vy(i,j,0);
683 
684  Real met_h_zeta;
685  met_h_zeta = Compute_h_zeta_AtEdgeCenterI(i,j,k,dxInv,z_nd);
686 
687  tau23(i,j,k) = 0.5 * ( (v(i, j, k) - v(i, j , k-1))*dxInv[2]/met_h_zeta
688  + (w(i, j, k) - w(i, j-1, k ))*dxInv[1]*mfy );
689  tau32(i,j,k) = tau23(i,j,k);
690  });
691  }
692 
693  //***********************************************************************************
694  // Fill the interior cells
695  //***********************************************************************************
696  // Cell centered strains
697  ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
698  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)
699  - 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) );
700  Real idz0 = 1.0 / dz0;
701 
702  Real GradUz = (k == 0) ?
703  idz0 * ( u(i ,j ,k+1) + u(i-1,j ,k+1)
704  - u(i ,j ,k ) - u(i-1,j ,k ) ) :
705  0.5 * idz0 * ( u(i ,j ,k+1) + u(i-1,j ,k+1)
706  - u(i ,j ,k-1) - u(i-1,j ,k-1) );
707  Real GradVz = (k == 0) ?
708  idz0 * ( v(i ,j ,k+1) + v(i ,j-1,k+1)
709  - v(i ,j ,k ) - v(i ,j-1,k ) ) :
710  0.5 * idz0 * ( v(i ,j ,k+1) + v(i ,j-1,k+1)
711  - v(i ,j ,k-1) - v(i ,j-1,k-1) );
712 
713  Real mfx = mf_mx(i,j,0);
714  Real mfy = mf_my(i,j,0);
715 
716  Real met_h_xi,met_h_eta,met_h_zeta;
717  met_h_xi = Compute_h_xi_AtCellCenter (i,j,k,dxInv,z_nd);
718  met_h_eta = Compute_h_eta_AtCellCenter (i,j,k,dxInv,z_nd);
719  met_h_zeta = detJ(i,j,k);
720 
721  tau11(i,j,k) = ( (u(i+1, j, k) - u(i, j, k))*dxInv[0] - met_h_xi*GradUz ) * mfx;
722  tau22(i,j,k) = ( (v(i, j+1, k) - v(i, j, k))*dxInv[1] - met_h_eta*GradVz ) * mfy;
723  tau33(i,j,k) = ( w(i, j, k+1) - w(i, j, k) )*dxInv[2]/met_h_zeta;
724  });
725 
726  // Off-diagonal strains
727  ParallelFor(tbxxy,tbxxz,tbxyz,
728  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
729  Real dz0 = ( z_nd(i,j,k+1) - z_nd(i,j,k-1) );
730  Real idz0 = 1.0 / dz0;
731 
732  Real GradUz = (k == 0) ?
733  idz0 * ( u(i ,j ,k+1) + u(i ,j-1,k+1)
734  - u(i ,j ,k ) - u(i ,j-1,k ) ) :
735  0.5 * idz0 * ( u(i ,j ,k+1) + u(i ,j-1,k+1)
736  - u(i ,j ,k-1) - u(i ,j-1,k-1) );
737  Real GradVz = (k == 0) ?
738  idz0 * ( v(i ,j ,k+1) + v(i-1,j ,k+1)
739  - v(i ,j ,k ) - v(i-1,j ,k ) ) :
740  0.5 * idz0 * ( v(i ,j ,k+1) + v(i-1,j ,k+1)
741  - v(i ,j ,k-1) - v(i-1,j ,k-1) );
742 
743  Real mfy = 0.5 * (mf_uy(i,j,0) + mf_uy(i ,j-1,0));
744  Real mfx = 0.5 * (mf_vx(i,j,0) + mf_vx(i-1,j ,0));
745 
746  Real met_h_xi,met_h_eta;
747  met_h_xi = Compute_h_xi_AtEdgeCenterK (i,j,k,dxInv,z_nd);
748  met_h_eta = Compute_h_eta_AtEdgeCenterK (i,j,k,dxInv,z_nd);
749 
750  tau12(i,j,k) = 0.5 * ( (u(i, j, k) - u(i , j-1, k))*dxInv[1]*mfy
751  + (v(i, j, k) - v(i-1, j , k))*dxInv[0]*mfx
752  - (met_h_eta)*GradUz*mfy
753  - (met_h_xi)*GradVz*mfx );
754  tau21(i,j,k) = tau12(i,j,k);
755  },
756  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
757  Real dz0 = 0.5 * ( z_nd(i,j,k+1) + z_nd(i,j+1,k+1)
758  - z_nd(i,j,k-1) - z_nd(i,j+1,k-1) );
759  Real idz0 = 1.0 / dz0;
760 
761  Real GradWz = 0.5 * idz0 * ( w(i ,j ,k+1) + w(i-1,j ,k+1)
762  - w(i ,j ,k-1) - w(i-1,j ,k-1) );
763 
764  Real mfx = mf_ux(i,j,0);
765 
766  Real met_h_xi,met_h_zeta;
767  met_h_xi = Compute_h_xi_AtEdgeCenterJ (i,j,k,dxInv,z_nd);
768  met_h_zeta = Compute_h_zeta_AtEdgeCenterJ(i,j,k,dxInv,z_nd);
769 
770  tau13(i,j,k) = 0.5 * ( (u(i, j, k) - u(i , j, k-1))*dxInv[2]/met_h_zeta
771  + ( (w(i, j, k) - w(i-1, j, k ))*dxInv[0]
772  - (met_h_xi)*GradWz ) * mfx );
773  tau31(i,j,k) = tau13(i,j,k);
774  },
775  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
776  Real dz0 = 0.5 * ( z_nd(i,j,k+1) + z_nd(i+1,j,k+1)
777  - z_nd(i,j,k-1) - z_nd(i+1,j,k-1) );
778  Real idz0 = 1.0 / dz0;
779 
780  Real GradWz = 0.5 * idz0 * ( w(i ,j ,k+1) + w(i ,j-1,k+1)
781  - w(i ,j ,k-1) - w(i ,j-1,k-1) );
782 
783  Real mfy = mf_vy(i,j,0);
784 
785  Real met_h_eta,met_h_zeta;
786  met_h_eta = Compute_h_eta_AtEdgeCenterI (i,j,k,dxInv,z_nd);
787  met_h_zeta = Compute_h_zeta_AtEdgeCenterI(i,j,k,dxInv,z_nd);
788 
789  tau23(i,j,k) = 0.5 * ( (v(i, j, k) - v(i, j , k-1))*dxInv[2]/met_h_zeta
790  + ( (w(i, j, k) - w(i, j-1, k ))*dxInv[1]
791  - (met_h_eta)*GradWz ) * mfy );
792  tau32(i,j,k) = tau23(i,j,k);
793  });
794 
795  if (SmnSmn_a) {
796  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
797  {
798  SmnSmn_a(i,j,k) = ComputeSmnSmn(i,j,k,
799  tau11,tau22,tau33,
800  tau12,tau13,tau23);
801  });
802  }
803 }
tau12
@ tau12
Definition: ERF_DataStruct.H:30
tau23
@ tau23
Definition: ERF_DataStruct.H:30
tau33
@ tau33
Definition: ERF_DataStruct.H:30
tau22
@ tau22
Definition: ERF_DataStruct.H:30
tau11
@ tau11
Definition: ERF_DataStruct.H:30
tau32
@ tau32
Definition: ERF_DataStruct.H:30
tau31
@ tau31
Definition: ERF_DataStruct.H:30
tau21
@ tau21
Definition: ERF_DataStruct.H:30
tau13
@ tau13
Definition: ERF_DataStruct.H:30
ComputeSmnSmn
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real ComputeSmnSmn(int &i, int &j, int &k, const amrex::Array4< amrex::Real const > &tau11, const amrex::Array4< amrex::Real const > &tau22, const amrex::Array4< amrex::Real const > &tau33, const amrex::Array4< amrex::Real const > &tau12, const amrex::Array4< amrex::Real const > &tau13, const amrex::Array4< amrex::Real const > &tau23)
Definition: ERF_EddyViscosity.H:63
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: