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

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: