ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
Functions
ERF_ComputeStrain_S.cpp File Reference
#include <ERF_Diffusion.H>
#include "ERF_EddyViscosity.H"
#include <ERF_TerrainMetrics.H>
Include dependency graph for ERF_ComputeStrain_S.cpp:

Functions

void ComputeStrain_S (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 Gpu::DeviceVector< Real > &stretched_dz_d, 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_S()

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