ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
Functions
ERF_ComputeStrain_T.cpp File Reference
#include <ERF_Diffusion.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)
 

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 
)

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
56 {
57  // Convert domain to each index type to test if we are on dirichlet boundary
58  Box domain_xy = convert(domain, tbxxy.ixType());
59  Box domain_xz = convert(domain, tbxxz.ixType());
60  Box domain_yz = convert(domain, tbxyz.ixType());
61 
62  const auto& dom_lo = lbound(domain);
63  const auto& dom_hi = ubound(domain);
64 
65  // Dirichlet on left or right plane
66  bool xl_v_dir = ( (bc_ptr[BCVars::yvel_bc].lo(0) == ERFBCType::ext_dir) ||
67  (bc_ptr[BCVars::yvel_bc].lo(0) == ERFBCType::ext_dir_upwind) ||
68  (bc_ptr[BCVars::yvel_bc].lo(0) == ERFBCType::ext_dir_ingested) );
69  xl_v_dir = ( xl_v_dir && (tbxxy.smallEnd(0) == domain_xy.smallEnd(0)) );
70 
71  bool xh_v_dir = ( (bc_ptr[BCVars::yvel_bc].hi(0) == ERFBCType::ext_dir) ||
72  (bc_ptr[BCVars::yvel_bc].hi(0) == ERFBCType::ext_dir_upwind) ||
73  (bc_ptr[BCVars::yvel_bc].hi(0) == ERFBCType::ext_dir_ingested) );
74  xh_v_dir = ( xh_v_dir && (tbxxy.bigEnd(0) == domain_xy.bigEnd(0)) );
75 
76  bool xl_w_dir = ( (bc_ptr[BCVars::zvel_bc].lo(0) == ERFBCType::ext_dir) ||
77  (bc_ptr[BCVars::zvel_bc].lo(0) == ERFBCType::ext_dir_upwind) ||
78  (bc_ptr[BCVars::zvel_bc].lo(0) == ERFBCType::ext_dir_ingested) );
79  xl_w_dir = ( xl_w_dir && (tbxxz.smallEnd(0) == domain_xz.smallEnd(0)) );
80 
81  bool xh_w_dir = ( (bc_ptr[BCVars::zvel_bc].hi(0) == ERFBCType::ext_dir) ||
82  (bc_ptr[BCVars::zvel_bc].hi(0) == ERFBCType::ext_dir_upwind) ||
83  (bc_ptr[BCVars::zvel_bc].hi(0) == ERFBCType::ext_dir_ingested) );
84  xh_w_dir = ( xh_w_dir && (tbxxz.bigEnd(0) == domain_xz.bigEnd(0)) );
85 
86  // Dirichlet on front or back plane
87  bool yl_u_dir = ( (bc_ptr[BCVars::xvel_bc].lo(1) == ERFBCType::ext_dir) ||
88  (bc_ptr[BCVars::xvel_bc].lo(1) == ERFBCType::ext_dir_upwind) ||
89  (bc_ptr[BCVars::xvel_bc].lo(1) == ERFBCType::ext_dir_ingested) );
90  yl_u_dir = ( yl_u_dir && (tbxxy.smallEnd(1) == domain_xy.smallEnd(1)) );
91 
92  bool yh_u_dir = ( (bc_ptr[BCVars::xvel_bc].hi(1) == ERFBCType::ext_dir) ||
93  (bc_ptr[BCVars::xvel_bc].hi(1) == ERFBCType::ext_dir_upwind) ||
94  (bc_ptr[BCVars::xvel_bc].hi(1) == ERFBCType::ext_dir_ingested) );
95  yh_u_dir = ( yh_u_dir && (tbxxy.bigEnd(1) == domain_xy.bigEnd(1)) );
96 
97  bool yl_w_dir = ( (bc_ptr[BCVars::zvel_bc].lo(1) == ERFBCType::ext_dir) ||
98  (bc_ptr[BCVars::zvel_bc].lo(1) == ERFBCType::ext_dir_upwind) ||
99  (bc_ptr[BCVars::zvel_bc].lo(1) == ERFBCType::ext_dir_ingested) );
100  yl_w_dir = ( yl_w_dir && (tbxyz.smallEnd(1) == domain_yz.smallEnd(1)) );
101 
102  bool yh_w_dir = ( (bc_ptr[BCVars::zvel_bc].hi(1) == ERFBCType::ext_dir) ||
103  (bc_ptr[BCVars::zvel_bc].hi(1) == ERFBCType::ext_dir_upwind) ||
104  (bc_ptr[BCVars::zvel_bc].hi(1) == ERFBCType::ext_dir_ingested) );
105  yh_w_dir = ( yh_w_dir && (tbxyz.bigEnd(1) == domain_yz.bigEnd(1)) );
106 
107  // Dirichlet on top or bottom plane
108  bool zl_u_dir = ( (bc_ptr[BCVars::xvel_bc].lo(2) == ERFBCType::ext_dir) ||
109  (bc_ptr[BCVars::xvel_bc].lo(2) == ERFBCType::ext_dir_ingested) );
110  zl_u_dir = ( zl_u_dir && (tbxxz.smallEnd(2) == domain_xz.smallEnd(2)) );
111 
112  bool zh_u_dir = ( (bc_ptr[BCVars::xvel_bc].hi(2) == ERFBCType::ext_dir) ||
113  (bc_ptr[BCVars::xvel_bc].hi(2) == ERFBCType::ext_dir_ingested) );
114  zh_u_dir = ( zh_u_dir && (tbxxz.bigEnd(2) == domain_xz.bigEnd(2)) );
115 
116  bool zl_v_dir = ( (bc_ptr[BCVars::yvel_bc].lo(2) == ERFBCType::ext_dir) ||
117  (bc_ptr[BCVars::yvel_bc].lo(2) == ERFBCType::ext_dir_ingested) );
118  zl_v_dir = ( zl_v_dir && (tbxyz.smallEnd(2) == domain_yz.smallEnd(2)) );
119 
120  bool zh_v_dir = ( (bc_ptr[BCVars::yvel_bc].hi(2) == ERFBCType::ext_dir) ||
121  (bc_ptr[BCVars::yvel_bc].hi(2) == ERFBCType::ext_dir_ingested) );
122  zh_v_dir = ( zh_v_dir && (tbxyz.bigEnd(2) == domain_yz.bigEnd(2)) );
123 
124  //***********************************************************************************
125  // X-Dirichlet
126  //***********************************************************************************
127  if (xl_v_dir)
128  {
129  Box planexy = tbxxy; planexy.setBig(0, planexy.smallEnd(0) );
130  tbxxy.growLo(0,-1);
131  bool need_to_test = (bc_ptr[BCVars::yvel_bc].lo(0) == ERFBCType::ext_dir_upwind) ? true : false;
132 
133  ParallelFor(planexy,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
134  Real inv_dist = (k == 0) ? Real(2.0) / (z_nd(i,j,k+2) - z_nd(i,j,k )) :
135  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));
136 
137  Real GradUz = (k == 0) ?
138  inv_dist * ( u(i ,j ,k+1) + u(i ,j-1,k+1)
139  - u(i ,j ,k ) - u(i ,j-1,k ) ) :
140  inv_dist * ( u(i ,j ,k+1) + u(i ,j-1,k+1)
141  - u(i ,j ,k-1) - u(i ,j-1,k-1) );
142  Real GradVz = (k == 0) ?
143  inv_dist * ( v(i ,j ,k+1) + v(i-1,j ,k+1)
144  - v(i ,j ,k ) - v(i-1,j ,k ) ) :
145  inv_dist * ( v(i ,j ,k+1) + v(i-1,j ,k+1)
146  - v(i ,j ,k-1) - v(i-1,j ,k-1) );
147 
148  Real mfy = 0.5 * (mf_uy(i,j,0) + mf_uy(i ,j-1,0));
149  Real mfx = 0.5 * (mf_vx(i,j,0) + mf_vx(i-1,j ,0));
150 
151  Real met_h_xi = Compute_h_xi_AtEdgeCenterK (i,j,k,dxInv,z_nd);
152  Real met_h_eta = Compute_h_eta_AtEdgeCenterK (i,j,k,dxInv,z_nd);
153 
154  if (!need_to_test || u(dom_lo.x,j,k) >= 0.) {
155  tau12(i,j,k) = 0.5 * ( (u(i, j, k) - u(i, j-1, k) )*dxInv[1]*mfy
156  + (-(8./3.) * v(i-1,j,k) + 3. * v(i,j,k) - (1./3.) * v(i+1,j,k))*dxInv[0]*mfx
157  - (met_h_eta)*GradUz*mfy
158  - (met_h_xi )*GradVz*mfx );
159  } else {
160  tau12(i,j,k) = 0.5 * ( (u(i, j, k) - u(i , j-1, k))*dxInv[1]*mfy
161  + (v(i, j, k) - v(i-1, j , k))*dxInv[0]*mfx
162  - (met_h_eta)*GradUz*mfy
163  - (met_h_xi )*GradVz*mfx );
164  }
165  tau21(i,j,k) = tau12(i,j,k);
166  });
167  }
168  if (xh_v_dir)
169  {
170  Box planexy = tbxxy; planexy.setSmall(0, planexy.bigEnd(0) );
171  tbxxy.growHi(0,-1);
172  bool need_to_test = (bc_ptr[BCVars::yvel_bc].hi(0) == ERFBCType::ext_dir_upwind) ? true : false;
173 
174  ParallelFor(planexy,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
175  Real inv_dist = (k == 0) ? Real(2.0) / (z_nd(i,j,k+2) - z_nd(i,j,k )) :
176  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));
177 
178  Real GradUz = (k == 0) ?
179  inv_dist * ( u(i ,j ,k+1) + u(i ,j-1,k+1)
180  - u(i ,j ,k ) - u(i ,j-1,k ) ) :
181  inv_dist * ( u(i ,j ,k+1) + u(i ,j-1,k+1)
182  - u(i ,j ,k-1) - u(i ,j-1,k-1) );
183  Real GradVz = (k == 0) ?
184  inv_dist * ( v(i ,j ,k+1) + v(i-1,j ,k+1)
185  - v(i ,j ,k ) - v(i-1,j ,k ) ) :
186  inv_dist * ( v(i ,j ,k+1) + v(i-1,j ,k+1)
187  - v(i ,j ,k-1) - v(i-1,j ,k-1) );
188 
189  Real mfy = 0.5 * (mf_uy(i,j,0) + mf_uy(i ,j-1,0));
190  Real mfx = 0.5 * (mf_vx(i,j,0) + mf_vx(i-1,j ,0));
191 
192  Real met_h_xi = Compute_h_xi_AtEdgeCenterK (i,j,k,dxInv,z_nd);
193  Real met_h_eta = Compute_h_eta_AtEdgeCenterK (i,j,k,dxInv,z_nd);
194 
195  if (!need_to_test || u(dom_hi.x+1,j,k) <= 0.) {
196  tau12(i,j,k) = 0.5 * ( (u(i, j, k) - u(i, j-1, k) )*dxInv[1]*mfy
197  - (-(8./3.) * v(i,j,k) + 3. * v(i-1,j,k) - (1./3.) * v(i-2,j,k))*dxInv[0]*mfx
198  - (met_h_eta)*GradUz*mfy
199  - (met_h_xi )*GradVz*mfx );
200  } else {
201  tau12(i,j,k) = 0.5 * ( (u(i, j, k) - u(i , j-1, k))*dxInv[1]*mfy
202  + (v(i, j, k) - v(i-1, j , k))*dxInv[0]*mfx
203  - (met_h_eta)*GradUz*mfy
204  - (met_h_xi )*GradVz*mfx );
205  }
206  tau21(i,j,k) = tau12(i,j,k);
207  });
208  }
209 
210  if (xl_w_dir)
211  {
212  Box planexz = tbxxz; planexz.setBig(0, planexz.smallEnd(0) );
213  planexz.setSmall(2, planexz.smallEnd(2)+1 ); planexz.setBig(2, planexz.bigEnd(2)-1 );
214  tbxxz.growLo(0,-1);
215  bool need_to_test = (bc_ptr[BCVars::zvel_bc].lo(0) == ERFBCType::ext_dir_upwind) ? true : false;
216 
217  ParallelFor(planexz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
218  Real dz0 = 0.5 * ( z_nd(i,j,k+1) + z_nd(i,j+1,k+1)
219  - z_nd(i,j,k-1) - z_nd(i,j+1,k-1) );
220  Real idz0 = 1.0 / dz0;
221 
222  Real GradWz = 0.5 * idz0 * ( w(i ,j ,k+1) + w(i-1,j ,k+1)
223  - w(i ,j ,k-1) - w(i-1,j ,k-1) );
224  Real mfx = mf_ux(i,j,0);
225 
226  Real met_h_xi = Compute_h_xi_AtEdgeCenterJ (i,j,k,dxInv,z_nd);
227  Real met_h_zeta = Compute_h_zeta_AtEdgeCenterJ(i,j,k,dxInv,z_nd);
228 
229  if (!need_to_test || u(dom_lo.x,j,k) <= 0.) {
230  tau13(i,j,k) = 0.5 * ( (u(i, j, k) - u(i, j, k-1))*dxInv[2]/met_h_zeta
231  + ( (-(8./3.) * w(i-1,j,k) + 3. * w(i,j,k) - (1./3.) * w(i+1,j,k))*dxInv[0]
232  - (met_h_xi)*GradWz ) * mfx );
233  } else {
234  tau13(i,j,k) = 0.5 * ( (u(i, j, k) - u(i , j, k-1))*dxInv[2]/met_h_zeta
235  + ( (w(i, j, k) - w(i-1, j, k ))*dxInv[0]
236  - (met_h_xi)*GradWz ) * mfx );
237  }
238  tau31(i,j,k) = tau13(i,j,k);
239  });
240  }
241 
242  if (xh_w_dir)
243  {
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  ParallelFor(planexz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
429  // Third order stencil with variable dz
430  Real dz0 = 0.5 * ( z_nd(i,j,k+1) + z_nd(i,j+1,k+1)
431  - z_nd(i,j,k ) - z_nd(i,j+1,k ) );
432  Real dz1 = 0.5 * ( z_nd(i,j,k+2) + z_nd(i,j+1,k+2)
433  - z_nd(i,j,k+1) - z_nd(i,j+1,k+1) );
434  Real idz0 = 1.0 / dz0;
435  Real f = (dz1 / dz0) + 2.0;
436  Real f2 = f*f;
437  Real c3 = 2.0 / (f - f2);
438  Real c2 = -f2*c3;
439  Real c1 = -(1.0-f2)*c3;
440 
441  Real GradWz = 0.5 * idz0 * ( w(i,j,k+1) + w(i-1,j,k+1)
442  - w(i,j,k ) - w(i-1,j,k ) );
443 
444  Real mfx = mf_ux(i,j,0);
445 
446  Real met_h_xi;
447  met_h_xi = Compute_h_xi_AtEdgeCenterJ (i,j,k,dxInv,z_nd);
448 
449  // Note: u(i,j,k) and u(i,j,k+1) are located at cell centers
450  // whereas u(i,j,k-1) is the value on the boundary
451  tau13(i,j,k) = 0.5 * ( (c1 * u(i,j,k-1) + c2 * u(i,j,k) + c3 * u(i,j,k+1))*idz0
452  + ( (w(i, j, k) - w(i-1, j, k))*dxInv[0]
453  - (met_h_xi)*GradWz ) * mfx );
454  tau31(i,j,k) = tau13(i,j,k);
455  });
456  }
457  if (zh_u_dir) {
458  // NOTE: h_xi = 0
459  Box planexz = tbxxz; planexz.setSmall(2, planexz.bigEnd(2) );
460  tbxxz.growHi(2,-1);
461  ParallelFor(planexz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
462  // Third order stencil with variable dz
463  Real dz0 = 0.5 * ( z_nd(i,j,k ) + z_nd(i,j+1,k )
464  - z_nd(i,j,k-1) - z_nd(i,j+1,k-1) );
465  Real dz1 = 0.5 * ( z_nd(i,j,k-1) + z_nd(i,j+1,k-1)
466  - z_nd(i,j,k-2) - z_nd(i,j+1,k-2) );
467  Real idz0 = 1.0 / dz0;
468  Real f = (dz1 / dz0) + 2.0;
469  Real f2 = f*f;
470  Real c3 = 2.0 / (f - f2);
471  Real c2 = -f2*c3;
472  Real c1 = -(1.0-f2)*c3;
473 
474  Real mfx = mf_ux(i,j,0);
475 
476  tau13(i,j,k) = 0.5 * ( -(c1 * u(i,j,k) + c2 * u(i,j,k-1) + c3 * u(i,j,k-2))*idz0
477  + (w(i, j, k) - w(i-1, j, k))*dxInv[0]*mfx );
478  tau31(i,j,k) = tau13(i,j,k);
479  });
480  }
481 
482  if (zl_v_dir) {
483  Box planeyz = tbxyz; planeyz.setBig(2, planeyz.smallEnd(2) );
484  tbxyz.growLo(2,-1);
485  ParallelFor(planeyz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
486  // Third order stencil with variable dz
487  Real dz0 = 0.5 * ( z_nd(i,j,k+1) + z_nd(i+1,j,k+1)
488  - z_nd(i,j,k ) - z_nd(i+1,j,k ) );
489  Real dz1 = 0.5 * ( z_nd(i,j,k+2) + z_nd(i+1,j,k+2)
490  - z_nd(i,j,k+1) - z_nd(i+1,j,k+1) );
491  Real idz0 = 1.0 / dz0;
492  Real f = (dz1 / dz0) + 2.0;
493  Real f2 = f*f;
494  Real c3 = 2.0 / (f - f2);
495  Real c2 = -f2*c3;
496  Real c1 = -(1.0-f2)*c3;
497 
498  Real GradWz = 0.5 * idz0 * ( w(i ,j ,k+1) + w(i ,j-1,k+1)
499  - w(i ,j ,k ) - w(i ,j-1,k ) );
500 
501  Real mfy = mf_vy(i,j,0);
502 
503  Real met_h_eta;
504  met_h_eta = Compute_h_eta_AtEdgeCenterI (i,j,k,dxInv,z_nd);
505 
506  // Note: v(i,j,k) and v(i,j,k+1) are located at cell centers
507  // whereas v(i,j,k-1) is the value on the boundary
508  tau23(i,j,k) = 0.5 * ( (c1 * v(i,j,k-1) + c2 * v(i,j,k ) + c3 * v(i,j,k+1))*idz0
509  + ( (w(i, j, k) - w(i, j-1, k))*dxInv[1]
510  - (met_h_eta)*GradWz ) * mfy );
511  tau32(i,j,k) = tau23(i,j,k);
512  });
513  }
514  if (zh_v_dir && (tbxyz.bigEnd(2) == domain_yz.bigEnd(2))) {
515  // NOTE: h_eta = 0
516  Box planeyz = tbxyz; planeyz.setSmall(2, planeyz.bigEnd(2) );
517  tbxyz.growHi(2,-1);
518  ParallelFor(planeyz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
519  // Third order stencil with variable dz
520  Real dz0 = 0.5 * ( z_nd(i,j,k ) + z_nd(i+1,j,k )
521  - z_nd(i,j,k-1) - z_nd(i+1,j,k-1) );
522  Real dz1 = 0.5 * ( z_nd(i,j,k-1) + z_nd(i+1,j,k-1)
523  - z_nd(i,j,k-2) - z_nd(i+1,j,k-2) );
524  Real idz0 = 1.0 / dz0;
525  Real f = (dz1 / dz0) + 2.0;
526  Real f2 = f*f;
527  Real c3 = 2.0 / (f - f2);
528  Real c2 = -f2*c3;
529  Real c1 = -(1.0-f2)*c3;
530 
531  Real mfy = mf_vy(i,j,0);
532 
533  tau23(i,j,k) = 0.5 * ( -(c1 * v(i,j,k ) + c2 * v(i,j,k-1) + c3 * v(i,j,k-2))*idz0
534  + (w(i, j, k) - w(i, j-1, k))*dxInv[1]*mfy );
535  tau32(i,j,k) = tau23(i,j,k);
536  });
537  }
538 
539  // HO derivatives w/ Dirichlet BC (\partival <var> / \partial z from terrain transform)
540  if (zl_u_dir && zl_v_dir) {
541  Box planecc = bxcc; planecc.setBig(2, planecc.smallEnd(2) );
542  bxcc.growLo(2,-1);
543  ParallelFor(planecc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
544  // Third order stencil with variable dz
545  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)
546  - z_nd(i,j,k ) - z_nd(i,j+1,k ) - z_nd(i+1,j,k ) - z_nd(i+1,j+1,k ) );
547  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)
548  - 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) );
549  Real idz0 = 1.0 / dz0;
550  Real f = (dz1 / dz0) + 2.0;
551  Real f2 = f*f;
552  Real c3 = 2.0 / (f - f2);
553  Real c2 = -f2*c3;
554  Real c1 = -(1.0-f2)*c3;
555 
556  Real GradUz = 0.5 * idz0 * ( (c1 * u(i ,j,k-1) + c2 * u(i ,j,k) + c3 * u(i ,j,k+1))
557  + (c1 * u(i-1,j,k-1) + c2 * u(i-1,j,k) + c3 * u(i-1,j,k+1)) );
558  Real GradVz = 0.5 * idz0 * ( (c1 * v(i,j ,k-1) + c2 * v(i,j ,k) + c3 * v(i,j ,k+1))
559  + (c1 * v(i,j-1,k-1) + c2 * v(i,j-1,k) + c3 * v(i,j-1,k+1)) );
560 
561  Real mfx = mf_mx(i,j,0);
562  Real mfy = mf_my(i,j,0);
563 
564  Real met_h_xi,met_h_eta;
565  met_h_xi = Compute_h_xi_AtCellCenter (i,j,k,dxInv,z_nd);
566  met_h_eta = Compute_h_eta_AtCellCenter (i,j,k,dxInv,z_nd);
567 
568  tau11(i,j,k) = ( (u(i+1, j, k) - u(i, j, k) )*dxInv[0]
569  - (met_h_xi)*GradUz ) * mfx;
570  tau22(i,j,k) = ( (v(i, j+1, k) - v(i, j, k) )*dxInv[1]
571  - (met_h_eta)*GradVz ) * mfy;
572  tau33(i,j,k) = ( w(i, j, k+1) - w(i, j, k) )*idz0;
573 
574  });
575 
576  Box planexy = tbxxy; planexy.setBig(2, planexy.smallEnd(2) );
577  tbxxy.growLo(2,-1);
578  ParallelFor(planexy,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
579  // Third order stencil with variable dz
580  Real dz0 = ( z_nd(i,j,k+1) - z_nd(i,j,k ) );
581  Real dz1 = ( z_nd(i,j,k+2) - z_nd(i,j,k+1) );
582  Real idz0 = 1.0 / dz0;
583  Real f = (dz1 / dz0) + 2.0;
584  Real f2 = f*f;
585  Real c3 = 2.0 / (f - f2);
586  Real c2 = -f2*c3;
587  Real c1 = -(1.0-f2)*c3;
588 
589  Real GradUz = 0.5 * idz0 * ( (c1 * u(i,j ,k-1) + c2 * u(i,j ,k) + c3 * u(i,j ,k+1))
590  + (c1 * u(i,j-1,k-1) + c2 * u(i,j-1,k) + c3 * u(i,j-1,k+1)) );
591  Real GradVz = 0.5 * idz0 * ( (c1 * v(i ,j,k-1) + c2 * v(i ,j,k) + c3 * v(i ,j,k+1))
592  + (c1 * v(i-1,j,k-1) + c2 * v(i-1,j,k) + c3 * v(i-1,j,k+1)) );
593 
594  Real mfy = 0.5 * (mf_uy(i,j,0) + mf_uy(i ,j-1,0));
595  Real mfx = 0.5 * (mf_vx(i,j,0) + mf_vx(i-1,j ,0));
596 
597  Real met_h_xi,met_h_eta;
598  met_h_xi = Compute_h_xi_AtEdgeCenterK (i,j,k,dxInv,z_nd);
599  met_h_eta = Compute_h_eta_AtEdgeCenterK (i,j,k,dxInv,z_nd);
600 
601  tau12(i,j,k) = 0.5 * ( (u(i, j, k) - u(i , j-1, k))*dxInv[1]*mfy
602  + (v(i, j, k) - v(i-1, j , k))*dxInv[0]*mfx
603  - (met_h_eta)*GradUz*mfy
604  - (met_h_xi )*GradVz*mfx );
605  tau21(i,j,k) = tau12(i,j,k);
606  });
607  }
608 
609  //***********************************************************************************
610  // Z-lo w/out Z-Dirichlet (GradWz extrapolation)
611  //***********************************************************************************
612  if (!zl_u_dir && (tbxxz.smallEnd(2) == domain_xz.smallEnd(2)) ) {
613  Box planexz = tbxxz; planexz.setBig(2, planexz.smallEnd(2));
614  tbxxz.growLo(2,-1);
615  ParallelFor(planexz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
616  Real GradWz = 0.5 * dxInv[2] * ( w(i ,j ,k+1) + w(i-1,j ,k+1)
617  - w(i ,j ,k ) - w(i-1,j ,k ) );
618 
619  Real mfx = mf_ux(i,j,0);
620 
621  Real met_h_xi,met_h_zeta;
622  met_h_xi = Compute_h_xi_AtEdgeCenterJ (i,j,k,dxInv,z_nd);
623  met_h_zeta = Compute_h_zeta_AtEdgeCenterJ(i,j,k,dxInv,z_nd);
624 
625  tau13(i,j,k) = 0.5 * ( (u(i, j, k) - u(i , j, k-1))*dxInv[2]/met_h_zeta
626  + ( (w(i, j, k) - w(i-1, j, k ))*dxInv[0]
627  - (met_h_xi/met_h_zeta)*GradWz ) * mfx);
628  tau31(i,j,k) = tau13(i,j,k);
629  });
630  }
631  if (!zl_v_dir && (tbxyz.smallEnd(2) == domain_yz.smallEnd(2))) {
632  Box planeyz = tbxyz; planeyz.setBig(2, planeyz.smallEnd(2) );
633  tbxyz.growLo(2,-1);
634  ParallelFor(planeyz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
635  Real GradWz = 0.5 * dxInv[2] * ( w(i ,j ,k+1) + w(i ,j-1,k+1)
636  - w(i ,j ,k ) - w(i ,j-1,k ) );
637 
638  Real mfy = mf_vy(i,j,0);
639 
640  Real met_h_eta,met_h_zeta;
641  met_h_eta = Compute_h_eta_AtEdgeCenterI (i,j,k,dxInv,z_nd);
642  met_h_zeta = Compute_h_zeta_AtEdgeCenterI(i,j,k,dxInv,z_nd);
643 
644  tau23(i,j,k) = 0.5 * ( (v(i, j, k) - v(i, j , k-1))*dxInv[2]/met_h_zeta
645  + ( (w(i, j, k) - w(i, j-1, k ))*dxInv[1]
646  - (met_h_eta/met_h_zeta)*GradWz ) * mfy );
647  tau32(i,j,k) = tau23(i,j,k);
648  });
649  }
650 
651  //***********************************************************************************
652  // Z-hi w/out Z-Dirichlet (h_xi = h_eta = 0)
653  //***********************************************************************************
654  if (!zh_u_dir && (tbxxz.bigEnd(2) == domain_xz.bigEnd(2))) {
655  Box planexz = tbxxz; planexz.setSmall(2, planexz.bigEnd(2) );
656  tbxxz.growHi(2,-1);
657  ParallelFor(planexz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
658  Real mfx = mf_ux(i,j,0);
659 
660  Real met_h_zeta;
661  met_h_zeta = Compute_h_zeta_AtEdgeCenterJ(i,j,k,dxInv,z_nd);
662 
663  tau13(i,j,k) = 0.5 * ( (u(i, j, k) - u(i , j, k-1))*dxInv[2]/met_h_zeta
664  + (w(i, j, k) - w(i-1, j, k ))*dxInv[0]*mfx );
665  tau31(i,j,k) = tau13(i,j,k);
666  });
667  }
668  if (!zh_v_dir && (tbxyz.bigEnd(2) == domain_yz.bigEnd(2))) {
669  Box planeyz = tbxyz; planeyz.setSmall(2, planeyz.bigEnd(2) );
670  tbxyz.growHi(2,-1);
671  ParallelFor(planeyz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
672  Real mfy = mf_vy(i,j,0);
673 
674  Real met_h_zeta;
675  met_h_zeta = Compute_h_zeta_AtEdgeCenterI(i,j,k,dxInv,z_nd);
676 
677  tau23(i,j,k) = 0.5 * ( (v(i, j, k) - v(i, j , k-1))*dxInv[2]/met_h_zeta
678  + (w(i, j, k) - w(i, j-1, k ))*dxInv[1]*mfy );
679  tau32(i,j,k) = tau23(i,j,k);
680  });
681  }
682 
683  //***********************************************************************************
684  // Fill the interior cells
685  //***********************************************************************************
686  // Cell centered strains
687  ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
688  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)
689  - 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) );
690  Real idz0 = 1.0 / dz0;
691 
692  Real GradUz = (k == 0) ?
693  idz0 * ( u(i ,j ,k+1) + u(i-1,j ,k+1)
694  - u(i ,j ,k ) - u(i-1,j ,k ) ) :
695  0.5 * idz0 * ( u(i ,j ,k+1) + u(i-1,j ,k+1)
696  - u(i ,j ,k-1) - u(i-1,j ,k-1) );
697  Real GradVz = (k == 0) ?
698  idz0 * ( v(i ,j ,k+1) + v(i ,j-1,k+1)
699  - v(i ,j ,k ) - v(i ,j-1,k ) ) :
700  0.5 * idz0 * ( v(i ,j ,k+1) + v(i ,j-1,k+1)
701  - v(i ,j ,k-1) - v(i ,j-1,k-1) );
702 
703  Real mfx = mf_mx(i,j,0);
704  Real mfy = mf_my(i,j,0);
705 
706  Real met_h_xi,met_h_eta,met_h_zeta;
707  met_h_xi = Compute_h_xi_AtCellCenter (i,j,k,dxInv,z_nd);
708  met_h_eta = Compute_h_eta_AtCellCenter (i,j,k,dxInv,z_nd);
709  met_h_zeta = detJ(i,j,k);
710 
711  tau11(i,j,k) = ( (u(i+1, j, k) - u(i, j, k))*dxInv[0] - met_h_xi*GradUz ) * mfx;
712  tau22(i,j,k) = ( (v(i, j+1, k) - v(i, j, k))*dxInv[1] - met_h_eta*GradVz ) * mfy;
713  tau33(i,j,k) = ( w(i, j, k+1) - w(i, j, k) )*dxInv[2]/met_h_zeta;
714  });
715 
716  // Off-diagonal strains
717  ParallelFor(tbxxy,tbxxz,tbxyz,
718  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
719  Real dz0 = ( z_nd(i,j,k+1) - z_nd(i,j,k-1) );
720  Real idz0 = 1.0 / dz0;
721 
722  Real GradUz = (k == 0) ?
723  idz0 * ( u(i ,j ,k+1) + u(i ,j-1,k+1)
724  - u(i ,j ,k ) - u(i ,j-1,k ) ) :
725  0.5 * idz0 * ( u(i ,j ,k+1) + u(i ,j-1,k+1)
726  - u(i ,j ,k-1) - u(i ,j-1,k-1) );
727  Real GradVz = (k == 0) ?
728  idz0 * ( v(i ,j ,k+1) + v(i-1,j ,k+1)
729  - v(i ,j ,k ) - v(i-1,j ,k ) ) :
730  0.5 * idz0 * ( v(i ,j ,k+1) + v(i-1,j ,k+1)
731  - v(i ,j ,k-1) - v(i-1,j ,k-1) );
732 
733  Real mfy = 0.5 * (mf_uy(i,j,0) + mf_uy(i ,j-1,0));
734  Real mfx = 0.5 * (mf_vx(i,j,0) + mf_vx(i-1,j ,0));
735 
736  Real met_h_xi,met_h_eta;
737  met_h_xi = Compute_h_xi_AtEdgeCenterK (i,j,k,dxInv,z_nd);
738  met_h_eta = Compute_h_eta_AtEdgeCenterK (i,j,k,dxInv,z_nd);
739 
740  tau12(i,j,k) = 0.5 * ( (u(i, j, k) - u(i , j-1, k))*dxInv[1]*mfy
741  + (v(i, j, k) - v(i-1, j , k))*dxInv[0]*mfx
742  - (met_h_eta)*GradUz*mfy
743  - (met_h_xi)*GradVz*mfx );
744  tau21(i,j,k) = tau12(i,j,k);
745  },
746  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
747  Real dz0 = 0.5 * ( z_nd(i,j,k+1) + z_nd(i,j+1,k+1)
748  - z_nd(i,j,k-1) - z_nd(i,j+1,k-1) );
749  Real idz0 = 1.0 / dz0;
750 
751  Real GradWz = 0.5 * idz0 * ( w(i ,j ,k+1) + w(i-1,j ,k+1)
752  - w(i ,j ,k-1) - w(i-1,j ,k-1) );
753 
754  Real mfx = mf_ux(i,j,0);
755 
756  Real met_h_xi,met_h_zeta;
757  met_h_xi = Compute_h_xi_AtEdgeCenterJ (i,j,k,dxInv,z_nd);
758  met_h_zeta = Compute_h_zeta_AtEdgeCenterJ(i,j,k,dxInv,z_nd);
759 
760  tau13(i,j,k) = 0.5 * ( (u(i, j, k) - u(i , j, k-1))*dxInv[2]/met_h_zeta
761  + ( (w(i, j, k) - w(i-1, j, k ))*dxInv[0]
762  - (met_h_xi)*GradWz ) * mfx );
763  tau31(i,j,k) = tau13(i,j,k);
764  },
765  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
766  Real dz0 = 0.5 * ( z_nd(i,j,k+1) + z_nd(i+1,j,k+1)
767  - z_nd(i,j,k-1) - z_nd(i+1,j,k-1) );
768  Real idz0 = 1.0 / dz0;
769 
770  Real GradWz = 0.5 * idz0 * ( w(i ,j ,k+1) + w(i ,j-1,k+1)
771  - w(i ,j ,k-1) - w(i ,j-1,k-1) );
772 
773  Real mfy = mf_vy(i,j,0);
774 
775  Real met_h_eta,met_h_zeta;
776  met_h_eta = Compute_h_eta_AtEdgeCenterI (i,j,k,dxInv,z_nd);
777  met_h_zeta = Compute_h_zeta_AtEdgeCenterI(i,j,k,dxInv,z_nd);
778 
779  tau23(i,j,k) = 0.5 * ( (v(i, j, k) - v(i, j , k-1))*dxInv[2]/met_h_zeta
780  + ( (w(i, j, k) - w(i, j-1, k ))*dxInv[1]
781  - (met_h_eta)*GradWz ) * mfy );
782  tau32(i,j,k) = tau23(i,j,k);
783  });
784 }
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
dom_hi
const auto & dom_hi
Definition: ERF_DiffSetup.H:10
dom_lo
const auto & dom_lo
Definition: ERF_DiffSetup.H:9
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:16
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:244
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:290
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:276
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:320
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:77
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:259
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:62
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:347
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:211

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: