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 > &tau13, Array4< Real > &tau21, Array4< Real > &tau23, Array4< Real > &tau31, Array4< Real > &tau32, const Array4< const Real > &z_nd, const Array4< const Real > &detJ, const BCRec *bc_ptr, const GpuArray< Real, AMREX_SPACEDIM > &dxInv, const Array4< const Real > &mf_m, const Array4< const Real > &mf_u, const Array4< const Real > &mf_v)
 

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 > &  tau13,
Array4< Real > &  tau21,
Array4< Real > &  tau23,
Array4< Real > &  tau31,
Array4< Real > &  tau32,
const Array4< const Real > &  z_nd,
const Array4< const Real > &  detJ,
const BCRec *  bc_ptr,
const GpuArray< Real, AMREX_SPACEDIM > &  dxInv,
const Array4< const Real > &  mf_m,
const Array4< const Real > &  mf_u,
const Array4< const Real > &  mf_v 
)

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_mmap factor at cell center
[in]mf_umap factor at x-face
[in]mf_vmap factor at y-face
45 {
46  // Convert domain to each index type to test if we are on dirichlet boundary
47  Box domain_xy = convert(domain, tbxxy.ixType());
48  Box domain_xz = convert(domain, tbxxz.ixType());
49  Box domain_yz = convert(domain, tbxyz.ixType());
50 
51  // Dirichlet on left or right plane
52  bool xl_v_dir = ( (bc_ptr[BCVars::yvel_bc].lo(0) == ERFBCType::ext_dir) ||
53  (bc_ptr[BCVars::yvel_bc].lo(0) == ERFBCType::ext_dir_ingested) );
54  xl_v_dir = ( xl_v_dir && (tbxxy.smallEnd(0) == domain_xy.smallEnd(0)) );
55 
56  bool xh_v_dir = ( (bc_ptr[BCVars::yvel_bc].hi(0) == ERFBCType::ext_dir) ||
57  (bc_ptr[BCVars::yvel_bc].hi(0) == ERFBCType::ext_dir_ingested) );
58  xh_v_dir = ( xh_v_dir && (tbxxy.bigEnd(0) == domain_xy.bigEnd(0)) );
59 
60  bool xl_w_dir = ( (bc_ptr[BCVars::zvel_bc].lo(0) == ERFBCType::ext_dir) ||
61  (bc_ptr[BCVars::zvel_bc].lo(0) == ERFBCType::ext_dir_ingested) );
62  xl_w_dir = ( xl_w_dir && (tbxxz.smallEnd(0) == domain_xz.smallEnd(0)) );
63 
64  bool xh_w_dir = ( (bc_ptr[BCVars::zvel_bc].hi(0) == ERFBCType::ext_dir) ||
65  (bc_ptr[BCVars::zvel_bc].hi(0) == ERFBCType::ext_dir_ingested) );
66  xh_w_dir = ( xh_w_dir && (tbxxz.bigEnd(0) == domain_xz.bigEnd(0)) );
67 
68  // Dirichlet on front or back plane
69  bool yl_u_dir = ( (bc_ptr[BCVars::xvel_bc].lo(1) == ERFBCType::ext_dir) ||
70  (bc_ptr[BCVars::xvel_bc].lo(1) == ERFBCType::ext_dir_ingested) );
71  yl_u_dir = ( yl_u_dir && (tbxxy.smallEnd(1) == domain_xy.smallEnd(1)) );
72 
73  bool yh_u_dir = ( (bc_ptr[BCVars::xvel_bc].hi(1) == ERFBCType::ext_dir) ||
74  (bc_ptr[BCVars::xvel_bc].hi(1) == ERFBCType::ext_dir_ingested) );
75  yh_u_dir = ( yh_u_dir && (tbxxy.bigEnd(1) == domain_xy.bigEnd(1)) );
76 
77  bool yl_w_dir = ( (bc_ptr[BCVars::zvel_bc].lo(1) == ERFBCType::ext_dir) ||
78  (bc_ptr[BCVars::zvel_bc].lo(1) == ERFBCType::ext_dir_ingested) );
79  yl_w_dir = ( yl_w_dir && (tbxyz.smallEnd(1) == domain_yz.smallEnd(1)) );
80 
81  bool yh_w_dir = ( (bc_ptr[BCVars::zvel_bc].hi(1) == ERFBCType::ext_dir) ||
82  (bc_ptr[BCVars::zvel_bc].hi(1) == ERFBCType::ext_dir_ingested) );
83  yh_w_dir = ( yh_w_dir && (tbxyz.bigEnd(1) == domain_yz.bigEnd(1)) );
84 
85  // Dirichlet on top or bottom plane
86  bool zl_u_dir = ( (bc_ptr[BCVars::xvel_bc].lo(2) == ERFBCType::ext_dir) ||
87  (bc_ptr[BCVars::xvel_bc].lo(2) == ERFBCType::ext_dir_ingested) );
88  zl_u_dir = ( zl_u_dir && (tbxxz.smallEnd(2) == domain_xz.smallEnd(2)) );
89 
90  bool zh_u_dir = ( (bc_ptr[BCVars::xvel_bc].hi(2) == ERFBCType::ext_dir) ||
91  (bc_ptr[BCVars::xvel_bc].hi(2) == ERFBCType::ext_dir_ingested) );
92  zh_u_dir = ( zh_u_dir && (tbxxz.bigEnd(2) == domain_xz.bigEnd(2)) );
93 
94  bool zl_v_dir = ( (bc_ptr[BCVars::yvel_bc].lo(2) == ERFBCType::ext_dir) ||
95  (bc_ptr[BCVars::yvel_bc].lo(2) == ERFBCType::ext_dir_ingested) );
96  zl_v_dir = ( zl_v_dir && (tbxyz.smallEnd(2) == domain_yz.smallEnd(2)) );
97 
98  bool zh_v_dir = ( (bc_ptr[BCVars::yvel_bc].hi(2) == ERFBCType::ext_dir) ||
99  (bc_ptr[BCVars::yvel_bc].hi(2) == ERFBCType::ext_dir_ingested) );
100  zh_v_dir = ( zh_v_dir && (tbxyz.bigEnd(2) == domain_yz.bigEnd(2)) );
101 
102  //***********************************************************************************
103  // X-Dirichlet
104  //***********************************************************************************
105  if (xl_v_dir) {
106  Box planexy = tbxxy; planexy.setBig(0, planexy.smallEnd(0) );
107  tbxxy.growLo(0,-1);
108  ParallelFor(planexy,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
109  Real GradUz = 0.25 * dxInv[2] * ( u(i ,j ,k+1)/mf_u(i,j,0) + u(i ,j-1,k+1)/mf_u(i,j-1,0)
110  -u(i ,j ,k-1)/mf_u(i,j,0) - u(i ,j-1,k-1)/mf_u(i,j-1,0) );
111  Real GradVz = 0.25 * dxInv[2] * ( v(i ,j ,k+1)/mf_v(i,j,0) + v(i-1,j ,k+1)/mf_v(i-1,j,0)
112  -v(i ,j ,k-1)/mf_v(i,j,0) - v(i-1,j ,k-1)/mf_v(i-1,j,0) );
113 
114  Real met_h_xi,met_h_eta,met_h_zeta;
115  met_h_xi = Compute_h_xi_AtEdgeCenterK (i,j,k,dxInv,z_nd);
116  met_h_eta = Compute_h_eta_AtEdgeCenterK (i,j,k,dxInv,z_nd);
117  met_h_zeta = Compute_h_zeta_AtEdgeCenterK(i,j,k,dxInv,z_nd);
118 
119  tau12(i,j,k) = 0.5 * ( (u(i, j, k)/mf_u(i,j,0) - u(i, j-1, k)/mf_u(i,j-1,0))*dxInv[1]
120  + (-(8./3.) * v(i-1,j,k)/mf_v(i-1,j,0) + 3. * v(i,j,k)/mf_v(i,j,0) - (1./3.) * v(i+1,j,k)/mf_v(i+1,j,0))*dxInv[0]
121  - (met_h_eta/met_h_zeta)*GradUz
122  - (met_h_xi /met_h_zeta)*GradVz ) * mf_u(i,j,0)*mf_u(i,j,0);
123  tau21(i,j,k) = tau12(i,j,k);
124  });
125  }
126  if (xh_v_dir) {
127  Box planexy = tbxxy; planexy.setSmall(0, planexy.bigEnd(0) );
128  tbxxy.growHi(0,-1);
129  ParallelFor(planexy,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
130  Real GradUz = 0.25 * dxInv[2] * ( u(i ,j ,k+1)/mf_u(i,j,0) + u(i ,j-1,k+1)/mf_u(i,j-1,0)
131  -u(i ,j ,k-1)/mf_u(i,j,0) - u(i ,j-1,k-1)/mf_u(i,j-1,0) );
132  Real GradVz = 0.25 * dxInv[2] * ( v(i ,j ,k+1)/mf_v(i,j,0) + v(i-1,j ,k+1)/mf_v(i-1,j,0)
133  -v(i ,j ,k-1)/mf_v(i,j,0) - v(i-1,j ,k-1)/mf_v(i-1,j,0) );
134 
135  Real met_h_xi,met_h_eta,met_h_zeta;
136  met_h_xi = Compute_h_xi_AtEdgeCenterK (i,j,k,dxInv,z_nd);
137  met_h_eta = Compute_h_eta_AtEdgeCenterK (i,j,k,dxInv,z_nd);
138  met_h_zeta = Compute_h_zeta_AtEdgeCenterK(i,j,k,dxInv,z_nd);
139 
140  tau12(i,j,k) = 0.5 * ( (u(i, j, k)/mf_u(i,j,0) - u(i, j-1, k)/mf_u(i,j-1,0))*dxInv[1]
141  - (-(8./3.) * v(i,j,k)/mf_v(i,j,0) + 3. * v(i-1,j,k)/mf_v(i-1,j,0) - (1./3.) * v(i-2,j,k)/mf_v(i-2,j,0))*dxInv[0]
142  - (met_h_eta/met_h_zeta)*GradUz
143  - (met_h_xi /met_h_zeta)*GradVz ) * mf_u(i,j,0)*mf_u(i,j,0);
144  tau21(i,j,k) = tau12(i,j,k);
145  });
146  }
147 
148  if (xl_w_dir) {
149  Box planexz = tbxxz; planexz.setBig(0, planexz.smallEnd(0) );
150  planexz.setSmall(2, planexz.smallEnd(2)+1 ); planexz.setBig(2, planexz.bigEnd(2)-1 );
151  tbxxz.growLo(0,-1);
152  ParallelFor(planexz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
153  Real GradWz = 0.25 * dxInv[2] * ( w(i ,j ,k+1) + w(i-1,j ,k+1)
154  - w(i ,j ,k-1) - w(i-1,j ,k-1) );
155 
156  Real met_h_xi,met_h_zeta;
157  met_h_xi = Compute_h_xi_AtEdgeCenterJ (i,j,k,dxInv,z_nd);
158  met_h_zeta = Compute_h_zeta_AtEdgeCenterJ(i,j,k,dxInv,z_nd);
159 
160  tau13(i,j,k) = 0.5 * ( (u(i, j, k) - u(i, j, k-1))*dxInv[2]/met_h_zeta
161  + ( (-(8./3.) * w(i-1,j,k) + 3. * w(i,j,k) - (1./3.) * w(i+1,j,k))*dxInv[0]
162  - (met_h_xi/met_h_zeta)*GradWz ) * mf_u(i,j,0) );
163  tau31(i,j,k) = tau13(i,j,k);
164  });
165  }
166  if (xh_w_dir) {
167  Box planexz = tbxxz; planexz.setSmall(0, planexz.bigEnd(0) );
168  planexz.setSmall(2, planexz.smallEnd(2)+1 ); planexz.setBig(2, planexz.bigEnd(2)-1 );
169  tbxxz.growHi(0,-1);
170  ParallelFor(planexz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
171  Real GradWz = 0.25 * dxInv[2] * ( w(i ,j ,k+1) + w(i-1,j ,k+1)
172  - w(i ,j ,k-1) - w(i-1,j ,k-1) );
173 
174  Real met_h_xi,met_h_zeta;
175  met_h_xi = Compute_h_xi_AtEdgeCenterJ (i,j,k,dxInv,z_nd);
176  met_h_zeta = Compute_h_zeta_AtEdgeCenterJ(i,j,k,dxInv,z_nd);
177 
178  tau13(i,j,k) = 0.5 * ( (u(i, j, k) - u(i, j, k-1))*dxInv[2]/met_h_zeta
179  - ( (-(8./3.) * w(i,j,k) + 3. * w(i-1,j,k) - (1./3.) * w(i-2,j,k))*dxInv[0]
180  - (met_h_xi/met_h_zeta)*GradWz ) * mf_u(i,j,0) );
181  tau31(i,j,k) = tau13(i,j,k);
182  });
183  }
184 
185  //***********************************************************************************
186  // Y-Dirichlet
187  //***********************************************************************************
188  if (yl_u_dir) {
189  Box planexy = tbxxy; planexy.setBig(1, planexy.smallEnd(1) );
190  tbxxy.growLo(1,-1);
191  ParallelFor(planexy,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
192  Real GradUz = 0.25 * dxInv[2] * ( u(i ,j ,k+1)/mf_u(i,j,0) + u(i ,j-1,k+1)/mf_u(i,j-1,0)
193  -u(i ,j ,k-1)/mf_u(i,j,0) - u(i ,j-1,k-1)/mf_u(i,j-1,0) );
194  Real GradVz = 0.25 * dxInv[2] * ( v(i ,j ,k+1)/mf_v(i,j,0) + v(i-1,j ,k+1)/mf_v(i-1,j,0)
195  -v(i ,j ,k-1)/mf_v(i,j,0) - v(i-1,j ,k-1)/mf_v(i-1,j,0) );
196 
197  Real met_h_xi,met_h_eta,met_h_zeta;
198  met_h_xi = Compute_h_xi_AtEdgeCenterK (i,j,k,dxInv,z_nd);
199  met_h_eta = Compute_h_eta_AtEdgeCenterK (i,j,k,dxInv,z_nd);
200  met_h_zeta = Compute_h_zeta_AtEdgeCenterK(i,j,k,dxInv,z_nd);
201 
202  tau12(i,j,k) = 0.5 * ( (-(8./3.) * u(i,j-1,k)/mf_u(i,j-1,0) + 3. * u(i,j,k)/mf_u(i,j,0) - (1./3.) * u(i,j+1,k)/mf_u(i,j+1,0))*dxInv[1]
203  + (v(i, j, k)/mf_v(i,j,0) - v(i-1, j, k)/mf_v(i,j,0))*dxInv[0]
204  - (met_h_eta/met_h_zeta)*GradUz
205  - (met_h_xi /met_h_zeta)*GradVz ) * mf_u(i,j,0)*mf_u(i,j,0);
206  tau21(i,j,k) = tau12(i,j,k);
207  });
208  }
209  if (yh_u_dir) {
210  Box planexy = tbxxy; planexy.setSmall(1, planexy.bigEnd(1) );
211  tbxxy.growHi(1,-1);
212  ParallelFor(planexy,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
213  Real GradUz = 0.25 * dxInv[2] * ( u(i ,j ,k+1)/mf_u(i,j,0) + u(i ,j-1,k+1)/mf_u(i,j-1,0)
214  -u(i ,j ,k-1)/mf_u(i,j,0) - u(i ,j-1,k-1)/mf_u(i,j-1,0) );
215  Real GradVz = 0.25 * dxInv[2] * ( v(i ,j ,k+1)/mf_v(i,j,0) + v(i-1,j ,k+1)/mf_v(i-1,j,0)
216  -v(i ,j ,k-1)/mf_v(i,j,0) - v(i-1,j ,k-1)/mf_v(i-1,j,0) );
217 
218  Real met_h_xi,met_h_eta,met_h_zeta;
219  met_h_xi = Compute_h_xi_AtEdgeCenterK (i,j,k,dxInv,z_nd);
220  met_h_eta = Compute_h_eta_AtEdgeCenterK (i,j,k,dxInv,z_nd);
221  met_h_zeta = Compute_h_zeta_AtEdgeCenterK(i,j,k,dxInv,z_nd);
222 
223  tau12(i,j,k) = 0.5 * ( -(-(8./3.) * u(i,j,k)/mf_u(i,j,0) + 3. * u(i,j-1,k)/mf_u(i,j-1,0) - (1./3.) * u(i,j-2,k)/mf_u(i,j-2,0))*dxInv[1] +
224  + (v(i, j, k)/mf_v(i,j,0) - v(i-1, j, k)/mf_v(i-1,j,0))*dxInv[0]
225  - (met_h_eta/met_h_zeta)*GradUz
226  - (met_h_xi /met_h_zeta)*GradVz ) * mf_u(i,j,0)*mf_u(i,j,0);
227  tau21(i,j,k) = tau12(i,j,k);
228  });
229  }
230 
231  if (yl_w_dir) {
232  Box planeyz = tbxyz; planeyz.setBig(1, planeyz.smallEnd(1) );
233  planeyz.setSmall(2, planeyz.smallEnd(2)+1 ); planeyz.setBig(2, planeyz.bigEnd(2)-1 );
234  tbxyz.growLo(1,-1);
235  ParallelFor(planeyz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
236  Real GradWz = 0.25 * dxInv[2] * ( w(i ,j ,k+1) + w(i ,j-1,k+1)
237  - w(i ,j ,k-1) - w(i ,j-1,k-1) );
238 
239  Real met_h_eta,met_h_zeta;
240  met_h_eta = Compute_h_eta_AtEdgeCenterI (i,j,k,dxInv,z_nd);
241  met_h_zeta = Compute_h_zeta_AtEdgeCenterI(i,j,k,dxInv,z_nd);
242 
243  tau23(i,j,k) = 0.5 * ( (v(i, j, k) - v(i, j, k-1))*dxInv[2]/met_h_zeta
244  + ( (-(8./3.) * w(i,j-1,k) + 3. * w(i,j ,k) - (1./3.) * w(i,j+1,k))*dxInv[1]*mf_v(i,j,0)
245  - (met_h_eta/met_h_zeta)*GradWz ) * mf_v(i,j,0) );
246  tau32(i,j,k) = tau23(i,j,k);
247  });
248  }
249  if (yh_w_dir) {
250  Box planeyz = tbxyz; planeyz.setSmall(1, planeyz.bigEnd(1) );
251  planeyz.setSmall(2, planeyz.smallEnd(2)+1 ); planeyz.setBig(2, planeyz.bigEnd(2)-1 );
252  tbxyz.growHi(1,-1);
253  ParallelFor(planeyz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
254  Real GradWz = 0.25 * dxInv[2] * ( w(i ,j ,k+1) + w(i ,j-1,k+1)
255  - w(i ,j ,k-1) - w(i ,j-1,k-1) );
256 
257  Real met_h_eta,met_h_zeta;
258  met_h_eta = Compute_h_eta_AtEdgeCenterI (i,j,k,dxInv,z_nd);
259  met_h_zeta = Compute_h_zeta_AtEdgeCenterI(i,j,k,dxInv,z_nd);
260 
261  tau23(i,j,k) = 0.5 * ( (v(i, j, k) - v(i, j, k-1))*dxInv[2]/met_h_zeta
262  - ( (-(8./3.) * w(i,j ,k) + 3. * w(i,j-1,k) - (1./3.) * w(i,j-2,k))*dxInv[1]
263  - (met_h_eta/met_h_zeta)*GradWz ) * mf_v(i,j,0) );
264  tau32(i,j,k) = tau23(i,j,k);
265  });
266  }
267 
268  //***********************************************************************************
269  // Z-Dirichlet
270  //***********************************************************************************
271  if (zl_u_dir) {
272  Box planexz = tbxxz; planexz.setBig(2, planexz.smallEnd(2) );
273  tbxxz.growLo(2,-1);
274  ParallelFor(planexz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
275  Real GradWz = 0.5 * dxInv[2] * ( w(i ,j ,k+1) + w(i-1,j ,k+1)
276  - w(i ,j ,k ) - w(i-1,j ,k ) );
277 
278  Real met_h_xi,met_h_zeta;
279  met_h_xi = Compute_h_xi_AtEdgeCenterJ (i,j,k,dxInv,z_nd);
280  met_h_zeta = Compute_h_zeta_AtEdgeCenterJ(i,j,k,dxInv,z_nd);
281 
282  tau13(i,j,k) = 0.5 * ( (-(8./3.) * u(i,j,k-1) + 3. * u(i,j,k) - (1./3.) * u(i,j,k+1))*dxInv[2]/met_h_zeta
283  + ( (w(i, j, k) - w(i-1, j, k))*dxInv[0]
284  - (met_h_xi/met_h_zeta)*GradWz ) * mf_u(i,j,0) );
285  tau31(i,j,k) = tau13(i,j,k);
286  });
287  }
288  if (zh_u_dir) {
289  Box planexz = tbxxz; planexz.setSmall(2, planexz.bigEnd(2) );
290  tbxxz.growHi(2,-1);
291  ParallelFor(planexz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
292  Real met_h_zeta;
293  met_h_zeta = Compute_h_zeta_AtEdgeCenterJ(i,j,k,dxInv,z_nd);
294 
295  tau13(i,j,k) = 0.5 * ( -(-(8./3.) * u(i,j,k) + 3. * u(i,j,k-1) - (1./3.) * u(i,j,k-2))*dxInv[2]/met_h_zeta
296  + (w(i, j, k) - w(i-1, j, k))*dxInv[0]*mf_u(i,j,0) );
297  tau31(i,j,k) = tau13(i,j,k);
298  });
299  }
300 
301  if (zl_v_dir) {
302  Box planeyz = tbxyz; planeyz.setBig(2, planeyz.smallEnd(2) );
303  tbxyz.growLo(2,-1);
304  ParallelFor(planeyz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
305  Real GradWz = 0.5 * dxInv[2] * ( w(i ,j ,k+1) + w(i ,j-1,k+1)
306  - w(i ,j ,k ) - w(i ,j-1,k ) );
307 
308  Real met_h_eta,met_h_zeta;
309  met_h_eta = Compute_h_eta_AtEdgeCenterI (i,j,k,dxInv,z_nd);
310  met_h_zeta = Compute_h_zeta_AtEdgeCenterI(i,j,k,dxInv,z_nd);
311 
312  tau23(i,j,k) = 0.5 * ( (-(8./3.) * v(i,j,k-1) + 3. * v(i,j,k ) - (1./3.) * v(i,j,k+1))*dxInv[2]/met_h_zeta
313  + ( (w(i, j, k) - w(i, j-1, k))*dxInv[1]
314  - (met_h_eta/met_h_zeta)*GradWz ) * mf_v(i,j,0) );
315  tau32(i,j,k) = tau23(i,j,k);
316  });
317  }
318  if (zh_v_dir && (tbxyz.bigEnd(2) == domain_yz.bigEnd(2))) {
319  Box planeyz = tbxyz; planeyz.setSmall(2, planeyz.bigEnd(2) );
320  tbxyz.growHi(2,-1);
321  ParallelFor(planeyz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
322  Real met_h_zeta;
323  met_h_zeta = Compute_h_zeta_AtEdgeCenterI(i,j,k,dxInv,z_nd);
324 
325  tau23(i,j,k) = 0.5 * ( -(-(8./3.) * v(i,j,k ) + 3. * v(i,j,k-1) - (1./3.) * v(i,j,k-2))*dxInv[2]/met_h_zeta
326  + (w(i, j, k) - w(i, j-1, k))*dxInv[1]*mf_v(i,j,0) );
327  tau32(i,j,k) = tau23(i,j,k);
328  });
329  }
330 
331  // HO derivatives w/ Dirichlet BC (\partival <var> / \partial z from terrain transform)
332  if (zl_u_dir && zl_v_dir) {
333  Box planecc = bxcc; planecc.setBig(2, planecc.smallEnd(2) );
334  bxcc.growLo(2,-1);
335  ParallelFor(planecc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
336  Real GradUz = 0.5 * dxInv[2] * ( (-(8./3.) * u(i ,j,k-1) + 3. * u(i ,j,k) - (1./3.) * u(i ,j,k+1))
337  + (-(8./3.) * u(i-1,j,k-1) + 3. * u(i-1,j,k) - (1./3.) * u(i-1,j,k+1)) );
338  Real GradVz = 0.5 * dxInv[2] * ( (-(8./3.) * v(i,j ,k-1) + 3. * v(i,j ,k) - (1./3.) * v(i,j ,k+1))
339  + (-(8./3.) * v(i,j-1,k-1) + 3. * v(i,j-1,k) - (1./3.) * v(i,j-1,k+1)) );
340 
341  Real met_h_xi,met_h_eta,met_h_zeta;
342  met_h_xi = Compute_h_xi_AtCellCenter (i,j,k,dxInv,z_nd);
343  met_h_eta = Compute_h_eta_AtCellCenter (i,j,k,dxInv,z_nd);
344  met_h_zeta = detJ(i,j,k);
345 
346  tau11(i,j,k) = ( (u(i+1, j, k)/mf_u(i+1,j,0) - u(i, j, k)/mf_u(i,j,0))*dxInv[0]
347  - (met_h_xi/met_h_zeta)*GradUz ) * mf_u(i,j,0)*mf_u(i,j,0);
348  tau22(i,j,k) = ( (v(i, j+1, k)/mf_v(i,j+1,0) - v(i, j, k)/mf_v(i,j,0))*dxInv[1]
349  - (met_h_eta/met_h_zeta)*GradVz ) * mf_v(i,j,0)*mf_v(i,j,0);
350  tau33(i,j,k) = (w(i, j, k+1) - w(i, j, k))*dxInv[2]/met_h_zeta;
351  });
352 
353  Box planexy = tbxxy; planexy.setBig(2, planexy.smallEnd(2) );
354  tbxxy.growLo(2,-1);
355  ParallelFor(planexy,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
356  Real GradUz = 0.5 * dxInv[2] * ( (-(8./3.) * u(i,j ,k-1) + 3. * u(i,j ,k) - (1./3.) * u(i,j ,k+1))
357  + (-(8./3.) * u(i,j-1,k-1) + 3. * u(i,j-1,k) - (1./3.) * u(i,j-1,k+1)) );
358  Real GradVz = 0.5 * dxInv[2] * ( (-(8./3.) * v(i ,j,k-1) + 3. * v(i ,j,k) - (1./3.) * v(i ,j,k+1))
359  + (-(8./3.) * v(i-1,j,k-1) + 3. * v(i-1,j,k) - (1./3.) * v(i-1,j,k+1)) );
360 
361  Real met_h_xi,met_h_eta,met_h_zeta;
362  met_h_xi = Compute_h_xi_AtEdgeCenterK (i,j,k,dxInv,z_nd);
363  met_h_eta = Compute_h_eta_AtEdgeCenterK (i,j,k,dxInv,z_nd);
364  met_h_zeta = Compute_h_zeta_AtEdgeCenterK(i,j,k,dxInv,z_nd);
365 
366  tau12(i,j,k) = 0.5 * ( (u(i, j, k)/mf_u(i,j,0) - u(i , j-1, k)/mf_u(i,j-1,0))*dxInv[1]
367  + (v(i, j, k)/mf_v(i,j,0) - v(i-1, j , k)/mf_v(i-1,j,0))*dxInv[0]
368  - (met_h_eta/met_h_zeta)*GradUz
369  - (met_h_xi /met_h_zeta)*GradVz ) * mf_u(i,j,0)*mf_u(i,j,0);
370  tau21(i,j,k) = tau12(i,j,k);
371  });
372  }
373 
374  //***********************************************************************************
375  // Z-lo w/out Z-Dirichlet (GradWz extrapolation)
376  //***********************************************************************************
377  if (!zl_u_dir && (tbxxz.smallEnd(2) == domain_xz.smallEnd(2)) ) {
378  Box planexz = tbxxz; planexz.setBig(2, planexz.smallEnd(2));
379  tbxxz.growLo(2,-1);
380  ParallelFor(planexz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
381  Real GradWz = 0.5 * dxInv[2] * ( w(i ,j ,k+1) + w(i-1,j ,k+1)
382  - w(i ,j ,k ) - w(i-1,j ,k ) );
383 
384  Real met_h_xi,met_h_zeta;
385  met_h_xi = Compute_h_xi_AtEdgeCenterJ (i,j,k,dxInv,z_nd);
386  met_h_zeta = Compute_h_zeta_AtEdgeCenterJ(i,j,k,dxInv,z_nd);
387 
388  tau13(i,j,k) = 0.5 * ( (u(i, j, k) - u(i , j, k-1))*dxInv[2]/met_h_zeta
389  + ( (w(i, j, k) - w(i-1, j, k ))*dxInv[0]
390  - (met_h_xi/met_h_zeta)*GradWz ) * mf_u(i,j,0) );
391  tau31(i,j,k) = tau13(i,j,k);
392  });
393  }
394  if (!zl_v_dir && (tbxyz.smallEnd(2) == domain_yz.smallEnd(2))) {
395  Box planeyz = tbxyz; planeyz.setBig(2, planeyz.smallEnd(2) );
396  tbxyz.growLo(2,-1);
397  ParallelFor(planeyz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
398  Real GradWz = 0.5 * dxInv[2] * ( w(i ,j ,k+1) + w(i ,j-1,k+1)
399  - w(i ,j ,k ) - w(i ,j-1,k ) );
400 
401  Real met_h_eta,met_h_zeta;
402  met_h_eta = Compute_h_eta_AtEdgeCenterI (i,j,k,dxInv,z_nd);
403  met_h_zeta = Compute_h_zeta_AtEdgeCenterI(i,j,k,dxInv,z_nd);
404 
405  tau23(i,j,k) = 0.5 * ( (v(i, j, k) - v(i, j , k-1))*dxInv[2]/met_h_zeta
406  + ( (w(i, j, k) - w(i, j-1, k ))*dxInv[1]
407  - (met_h_eta/met_h_zeta)*GradWz ) * mf_v(i,j,0) );
408  tau32(i,j,k) = tau23(i,j,k);
409  });
410  }
411 
412  //***********************************************************************************
413  // Z-hi w/out Z-Dirichlet (h_xi = h_eta = 0)
414  //***********************************************************************************
415  if (!zh_u_dir && (tbxxz.bigEnd(2) == domain_xz.bigEnd(2))) {
416  Box planexz = tbxxz; planexz.setSmall(2, planexz.bigEnd(2) );
417  tbxxz.growHi(2,-1);
418  ParallelFor(planexz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
419  Real met_h_zeta;
420  met_h_zeta = Compute_h_zeta_AtEdgeCenterJ(i,j,k,dxInv,z_nd);
421 
422  tau13(i,j,k) = 0.5 * ( (u(i, j, k) - u(i , j, k-1))*dxInv[2]/met_h_zeta
423  + (w(i, j, k) - w(i-1, j, k ))*dxInv[0]*mf_u(i,j,0) );
424  tau31(i,j,k) = tau13(i,j,k);
425  });
426  }
427  if (!zh_v_dir && (tbxyz.bigEnd(2) == domain_yz.bigEnd(2))) {
428  Box planeyz = tbxyz; planeyz.setSmall(2, planeyz.bigEnd(2) );
429  tbxyz.growHi(2,-1);
430  ParallelFor(planeyz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
431  Real met_h_zeta;
432  met_h_zeta = Compute_h_zeta_AtEdgeCenterI(i,j,k,dxInv,z_nd);
433 
434  tau23(i,j,k) = 0.5 * ( (v(i, j, k) - v(i, j , k-1))*dxInv[2]/met_h_zeta
435  + (w(i, j, k) - w(i, j-1, k ))*dxInv[1]*mf_v(i,j,0) );
436  tau32(i,j,k) = tau23(i,j,k);
437  });
438  }
439 
440  //***********************************************************************************
441  // Fill the interior cells
442  //***********************************************************************************
443  // Cell centered strains
444  ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
445  Real GradUz = 0.25 * dxInv[2] * ( u(i ,j ,k+1)/mf_u(i,j,0) + u(i-1,j ,k+1)/mf_u(i-1,j,0)
446  -u(i ,j ,k-1)/mf_u(i,j,0) - u(i-1,j ,k-1)/mf_u(i-1,j,0) );
447  Real GradVz = 0.25 * dxInv[2] * ( v(i ,j ,k+1)/mf_v(i,j,0) + v(i ,j-1,k+1)/mf_v(i,j-1,0)
448  -v(i ,j ,k-1)/mf_v(i,j,0) - v(i ,j-1,k-1)/mf_v(i,j-1,0) );
449 
450  Real met_h_xi,met_h_eta,met_h_zeta;
451  met_h_xi = Compute_h_xi_AtCellCenter (i,j,k,dxInv,z_nd);
452  met_h_eta = Compute_h_eta_AtCellCenter (i,j,k,dxInv,z_nd);
453  met_h_zeta = detJ(i,j,k);
454 
455  tau11(i,j,k) = ( (u(i+1, j, k)/mf_u(i+1,j,0) - u(i, j, k)/mf_u(i,j,0))*dxInv[0]
456  - (met_h_xi/met_h_zeta)*GradUz ) * mf_m(i,j,0)*mf_m(i,j,0);
457  tau22(i,j,k) = ( (v(i, j+1, k)/mf_v(i,j+1,0) - v(i, j, k)/mf_v(i,j,0))*dxInv[1]
458  - (met_h_eta/met_h_zeta)*GradVz ) * mf_m(i,j,0)*mf_m(i,j,0);
459  tau33(i,j,k) = (w(i, j, k+1) - w(i, j, k))*dxInv[2]/met_h_zeta;
460  });
461 
462  // Off-diagonal strains
463  ParallelFor(tbxxy,tbxxz,tbxyz,
464  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
465  Real GradUz = 0.25 * dxInv[2] * ( u(i ,j ,k+1)/mf_u(i,j,0) + u(i ,j-1,k+1)/mf_u(i,j-1,0)
466  -u(i ,j ,k-1)/mf_u(i,j,0) - u(i ,j-1,k-1)/mf_u(i,j-1,0) );
467  Real GradVz = 0.25 * dxInv[2] * ( v(i ,j ,k+1)/mf_v(i,j,0) + v(i-1,j ,k+1)/mf_v(i-1,j,0)
468  -v(i ,j ,k-1)/mf_v(i,j,0) - v(i-1,j ,k-1)/mf_v(i-1,j,0) );
469 
470  Real met_h_xi,met_h_eta,met_h_zeta;
471  met_h_xi = Compute_h_xi_AtEdgeCenterK (i,j,k,dxInv,z_nd);
472  met_h_eta = Compute_h_eta_AtEdgeCenterK (i,j,k,dxInv,z_nd);
473  met_h_zeta = Compute_h_zeta_AtEdgeCenterK(i,j,k,dxInv,z_nd);
474 
475  tau12(i,j,k) = 0.5 * ( (u(i, j, k)/mf_u(i,j,0) - u(i , j-1, k)/mf_u(i,j-1,0))*dxInv[1]
476  + (v(i, j, k)/mf_v(i,j,0) - v(i-1, j , k)/mf_v(i-1,j,0))*dxInv[0]
477  - (met_h_eta/met_h_zeta)*GradUz
478  - (met_h_xi /met_h_zeta)*GradVz ) * mf_u(i,j,0)*mf_u(i,j,0);
479  tau21(i,j,k) = tau12(i,j,k);
480  },
481  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
482  Real GradWz = 0.25 * dxInv[2] * ( w(i ,j ,k+1) + w(i-1,j ,k+1)
483  - w(i ,j ,k-1) - w(i-1,j ,k-1) );
484 
485  Real met_h_xi,met_h_zeta;
486  met_h_xi = Compute_h_xi_AtEdgeCenterJ (i,j,k,dxInv,z_nd);
487  met_h_zeta = Compute_h_zeta_AtEdgeCenterJ(i,j,k,dxInv,z_nd);
488 
489  tau13(i,j,k) = 0.5 * ( (u(i, j, k) - u(i , j, k-1))*dxInv[2]/met_h_zeta
490  + ( (w(i, j, k) - w(i-1, j, k ))*dxInv[0]
491  - (met_h_xi/met_h_zeta)*GradWz ) * mf_u(i,j,0) );
492  tau31(i,j,k) = tau13(i,j,k);
493  },
494  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
495  Real GradWz = 0.25 * dxInv[2] * ( w(i ,j ,k+1) + w(i ,j-1,k+1)
496  - w(i ,j ,k-1) - w(i ,j-1,k-1) );
497 
498  Real met_h_eta,met_h_zeta;
499  met_h_eta = Compute_h_eta_AtEdgeCenterI (i,j,k,dxInv,z_nd);
500  met_h_zeta = Compute_h_zeta_AtEdgeCenterI(i,j,k,dxInv,z_nd);
501 
502  tau23(i,j,k) = 0.5 * ( (v(i, j, k) - v(i, j , k-1))*dxInv[2]/met_h_zeta
503  + ( (w(i, j, k) - w(i, j-1, k ))*dxInv[1]
504  - (met_h_eta/met_h_zeta)*GradWz ) * mf_v(i,j,0) );
505  tau32(i,j,k) = tau23(i,j,k);
506  });
507 }
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:234
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:280
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:266
Compute_h_zeta_AtEdgeCenterK
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_zeta_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:221
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:310
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:69
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:249
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:54
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:337
BCVars::zvel_bc
@ zvel_bc
Definition: ERF_IndexDefines.H:90
BCVars::yvel_bc
@ yvel_bc
Definition: ERF_IndexDefines.H:89
BCVars::xvel_bc
@ xvel_bc
Definition: ERF_IndexDefines.H:88
ERFBCType::ext_dir_ingested
@ ext_dir_ingested
Definition: ERF_IndexDefines.H:183
ERFBCType::ext_dir
@ ext_dir
Definition: ERF_IndexDefines.H:180

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: