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 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 > &tau13i, Array4< Real > &tau23i)
 

Function Documentation

◆ ComputeStrain_S()

void ComputeStrain_S ( 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 > &  tau13i,
Array4< Real > &  tau23i 
)

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]tau13icontribution to strain from du/dz
[in]tau23icontribution to strain from dv/dz
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  Real du_dz = (u(i, j, k) - u(i, j, k-1))*dz_inv;
183  if (!need_to_test || u(dom_lo.x,j,k) <= 0.) {
184  tau13(i,j,k) = 0.5 * ( du_dz
185  + (-(8./3.) * w(i-1,j,k) + 3. * w(i,j,k) - (1./3.) * w(i+1,j,k))*dxInv[0] * mfx );
186  } else {
187  tau13(i,j,k) = 0.5 * ( du_dz
188  + (w(i, j, k) - w(i-1, j, k ))*dxInv[0] * mfx );
189  }
190  tau31(i,j,k) = tau13(i,j,k);
191 
192  if (tau13i) tau13i(i,j,k) = 0.5 * du_dz;
193  });
194  }
195 
196  if (xh_w_dir) {
197  Box planexz = tbxxz; planexz.setSmall(0, planexz.bigEnd(0) );
198  planexz.setSmall(2, planexz.smallEnd(2)+1 ); planexz.setBig(2, planexz.bigEnd(2)-1 );
199  tbxxz.growHi(0,-1);
200  bool need_to_test = (bc_ptr[BCVars::zvel_bc].hi(0) == ERFBCType::ext_dir_upwind) ? true : false;
201 
202  ParallelFor(planexz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
203  Real mfx = mf_ux(i,j,0);
204 
205  Real dz_inv = (k == 0) ? 1.0 / dz_ptr[k] : 2.0 / (dz_ptr[k] + dz_ptr[k-1]);
206 
207  Real du_dz = (u(i, j, k) - u(i, j, k-1))*dz_inv;
208  if (!need_to_test || u(dom_hi.x+1,j,k) <= 0.) {
209  tau13(i,j,k) = 0.5 * ( du_dz
210  - (-(8./3.) * w(i,j,k) + 3. * w(i-1,j,k) - (1./3.) * w(i-2,j,k))*dxInv[0] * mfx );
211  } else {
212  tau13(i,j,k) = 0.5 * ( du_dz
213  + (w(i, j, k) - w(i-1, j, k ))*dxInv[0] * mfx );
214  }
215  tau31(i,j,k) = tau13(i,j,k);
216 
217  if (tau13i) tau13i(i,j,k) = 0.5 * du_dz;
218  });
219  }
220 
221  //***********************************************************************************
222  // Y-Dirichlet
223  //***********************************************************************************
224  if (yl_u_dir) {
225  Box planexy = tbxxy; planexy.setBig(1, planexy.smallEnd(1) );
226  tbxxy.growLo(1,-1);
227  bool need_to_test = (bc_ptr[BCVars::xvel_bc].lo(1) == ERFBCType::ext_dir_upwind) ? true : false;
228 
229  ParallelFor(planexy,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
230  Real mfy = 0.5 * (mf_uy(i,j,0) + mf_uy(i ,j-1,0));
231  Real mfx = 0.5 * (mf_vx(i,j,0) + mf_vx(i-1,j ,0));
232 
233  if (!need_to_test || v(i,dom_lo.y,k) >= 0.) {
234  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
235  + (v(i, j, k) - v(i-1, j, k))*dxInv[0]*mfx);
236  } else {
237  tau12(i,j,k) = 0.5 * ( (u(i, j, k) - u(i , j-1, k))*dxInv[1]*mfy
238  + (v(i, j, k) - v(i-1, j , k))*dxInv[0]*mfx);
239  }
240  tau21(i,j,k) = tau12(i,j,k);
241  });
242  }
243  if (yh_u_dir) {
244  Box planexy = tbxxy; planexy.setSmall(1, planexy.bigEnd(1) );
245  tbxxy.growHi(1,-1);
246  bool need_to_test = (bc_ptr[BCVars::xvel_bc].hi(1) == ERFBCType::ext_dir_upwind) ? true : false;
247 
248  ParallelFor(planexy,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
249  Real mfy = 0.5 * (mf_uy(i,j,0) + mf_uy(i ,j-1,0));
250  Real mfx = 0.5 * (mf_vx(i,j,0) + mf_vx(i-1,j ,0));
251 
252  if (!need_to_test || v(i,dom_hi.y+1,k) >= 0.) {
253  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 +
254  + (v(i, j, k) - v(i-1, j, k))*dxInv[0]*mfx);
255  } else {
256  tau12(i,j,k) = 0.5 * ( (u(i, j, k) - u(i , j-1, k))*dxInv[1]*mfy
257  + (v(i, j, k) - v(i-1, j , k))*dxInv[0]*mfx);
258  }
259  tau21(i,j,k) = tau12(i,j,k);
260  });
261  }
262 
263  if (yl_w_dir) {
264  Box planeyz = tbxyz; planeyz.setBig(1, planeyz.smallEnd(1) );
265  planeyz.setSmall(2, planeyz.smallEnd(2)+1 ); planeyz.setBig(2, planeyz.bigEnd(2)-1 );
266  tbxyz.growLo(1,-1);
267  bool need_to_test = (bc_ptr[BCVars::zvel_bc].lo(1) == ERFBCType::ext_dir_upwind) ? true : false;
268 
269  ParallelFor(planeyz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
270  Real mfy = mf_vy(i,j,0);
271 
272  Real dz_inv = (k == 0) ? 1.0 / dz_ptr[k] : 2.0 / (dz_ptr[k] + dz_ptr[k-1]);
273 
274  Real dv_dz = (v(i, j, k) - v(i, j, k-1))*dz_inv;
275  if (!need_to_test || v(i,dom_lo.y,k) >= 0.) {
276  tau23(i,j,k) = 0.5 * ( dv_dz
277  + (-(8./3.) * w(i,j-1,k) + 3. * w(i,j ,k) - (1./3.) * w(i,j+1,k))*dxInv[1] * mfy );
278  } else {
279  tau23(i,j,k) = 0.5 * ( dv_dz
280  + (w(i, j, k) - w(i, j-1, k ))*dxInv[1] * mfy );
281  }
282  tau32(i,j,k) = tau23(i,j,k);
283 
284  if (tau23i) tau23i(i,j,k) = 0.5 * dv_dz;
285  });
286  }
287  if (yh_w_dir) {
288  Box planeyz = tbxyz; planeyz.setSmall(1, planeyz.bigEnd(1) );
289  planeyz.setSmall(2, planeyz.smallEnd(2)+1 ); planeyz.setBig(2, planeyz.bigEnd(2)-1 );
290  tbxyz.growHi(1,-1);
291  bool need_to_test = (bc_ptr[BCVars::zvel_bc].hi(1) == ERFBCType::ext_dir_upwind) ? true : false;
292 
293  ParallelFor(planeyz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
294  Real mfy = mf_vy(i,j,0);
295 
296  Real dz_inv = (k == 0) ? 1.0 / dz_ptr[k] : 2.0 / (dz_ptr[k] + dz_ptr[k-1]);
297 
298  Real dv_dz = (v(i, j, k) - v(i, j, k-1))*dz_inv;
299  if (!need_to_test || v(i,dom_hi.y+1,k) >= 0.) {
300  tau23(i,j,k) = 0.5 * ( dv_dz
301  - (-(8./3.) * w(i,j ,k) + 3. * w(i,j-1,k) - (1./3.) * w(i,j-2,k))*dxInv[1] * mfy );
302  } else {
303  tau23(i,j,k) = 0.5 * ( dv_dz
304  + (w(i, j, k) - w(i, j-1, k ))*dxInv[1] * mfy );
305  }
306  tau32(i,j,k) = tau23(i,j,k);
307 
308  if (tau23i) tau23i(i,j,k) = 0.5 * dv_dz;
309  });
310  }
311 
312  //***********************************************************************************
313  // Z-Dirichlet
314  //***********************************************************************************
315  if (zl_u_dir) {
316  Box planexz = tbxxz; planexz.setBig(2, planexz.smallEnd(2) );
317  tbxxz.growLo(2,-1);
318 
319  ParallelFor(planexz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
320  // Third order stencil with variable dz
321  Real dz0 = dz_ptr[k];
322  Real dz1 = dz_ptr[k+1];
323  Real idz0 = 1.0 / dz0;
324  Real f = (dz1 / dz0) + 2.0;
325  Real f2 = f*f;
326  Real c3 = 2.0 / (f - f2);
327  Real c2 = -f2*c3;
328  Real c1 = -(1.0-f2)*c3;
329 
330  Real mfx = mf_ux(i,j,0);
331 
332  Real du_dz = (c1 * u(i,j,k-1) + c2 * u(i,j,k) + c3 * u(i,j,k+1))*idz0;
333  tau13(i,j,k) = 0.5 * ( du_dz
334  + (w(i, j, k) - w(i-1, j, k))*dxInv[0] * mfx );
335  tau31(i,j,k) = tau13(i,j,k);
336 
337  if (tau13i) tau13i(i,j,k) = 0.5 * du_dz;
338  });
339  }
340  if (zh_u_dir) {
341  // NOTE: h_xi = 0
342  Box planexz = tbxxz; planexz.setSmall(2, planexz.bigEnd(2) );
343  tbxxz.growHi(2,-1);
344 
345  ParallelFor(planexz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
346  // Third order stencil with variable dz
347  Real dz0 = dz_ptr[k-1];
348  Real dz1 = dz_ptr[k-2];
349  Real idz0 = 1.0 / dz0;
350  Real f = (dz1 / dz0) + 2.0;
351  Real f2 = f*f;
352  Real c3 = 2.0 / (f - f2);
353  Real c2 = -f2*c3;
354  Real c1 = -(1.0-f2)*c3;
355 
356  Real mfx = mf_ux(i,j,0);
357 
358  Real du_dz = -(c1 * u(i,j,k) + c2 * u(i,j,k-1) + c3 * u(i,j,k-2))*idz0;
359  tau13(i,j,k) = 0.5 * ( du_dz
360  + (w(i, j, k) - w(i-1, j, k))*dxInv[0]*mfx );
361  tau31(i,j,k) = tau13(i,j,k);
362 
363  if (tau13i) tau13i(i,j,k) = 0.5 * du_dz;
364  });
365  }
366 
367  if (zl_v_dir) {
368  Box planeyz = tbxyz; planeyz.setBig(2, planeyz.smallEnd(2) );
369  tbxyz.growLo(2,-1);
370 
371  ParallelFor(planeyz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
372  // Third order stencil with variable dz
373  Real dz0 = dz_ptr[k];
374  Real dz1 = dz_ptr[k+1];
375  Real idz0 = 1.0 / dz0;
376  Real f = (dz1 / dz0) + 2.0;
377  Real f2 = f*f;
378  Real c3 = 2.0 / (f - f2);
379  Real c2 = -f2*c3;
380  Real c1 = -(1.0-f2)*c3;
381 
382  Real mfy = mf_vy(i,j,0);
383 
384  Real dv_dz = (c1 * v(i,j,k-1) + c2 * v(i,j,k ) + c3 * v(i,j,k+1))*idz0;
385  tau23(i,j,k) = 0.5 * ( dv_dz
386  + (w(i, j, k) - w(i, j-1, k))*dxInv[1] * mfy );
387  tau32(i,j,k) = tau23(i,j,k);
388 
389  if (tau23i) tau23i(i,j,k) = 0.5 * dv_dz;
390  });
391  }
392  if (zh_v_dir && (tbxyz.bigEnd(2) == domain_yz.bigEnd(2))) {
393  // NOTE: h_eta = 0
394  Box planeyz = tbxyz; planeyz.setSmall(2, planeyz.bigEnd(2) );
395  tbxyz.growHi(2,-1);
396 
397  ParallelFor(planeyz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
398  // Third order stencil with variable dz
399  Real dz0 = dz_ptr[k-1];
400  Real dz1 = dz_ptr[k-2];
401  Real idz0 = 1.0 / dz0;
402  Real f = (dz1 / dz0) + 2.0;
403  Real f2 = f*f;
404  Real c3 = 2.0 / (f - f2);
405  Real c2 = -f2*c3;
406  Real c1 = -(1.0-f2)*c3;
407 
408  Real mfy = mf_vy(i,j,0);
409 
410  Real dv_dz = -(c1 * v(i,j,k ) + c2 * v(i,j,k-1) + c3 * v(i,j,k-2))*idz0;
411  tau23(i,j,k) = 0.5 * ( dv_dz
412  + (w(i, j, k) - w(i, j-1, k))*dxInv[1]*mfy );
413  tau32(i,j,k) = tau23(i,j,k);
414 
415  if (tau23i) tau23i(i,j,k) = 0.5 * dv_dz;
416  });
417  }
418 
419  // HO derivatives w/ Dirichlet BC (\partival <var> / \partial z from terrain transform)
420  if (zl_u_dir && zl_v_dir) {
421  Box planecc = bxcc; planecc.setBig(2, planecc.smallEnd(2) );
422  bxcc.growLo(2,-1);
423 
424  ParallelFor(planecc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
425  Real dz0 = dz_ptr[k];
426  Real idz0 = 1.0 / dz0;
427 
428  Real mfx = mf_mx(i,j,0);
429  Real mfy = mf_my(i,j,0);
430 
431  tau11(i,j,k) = ( (u(i+1, j, k) - u(i, j, k) )*dxInv[0]) * mfx;
432  tau22(i,j,k) = ( (v(i, j+1, k) - v(i, j, k) )*dxInv[1]) * mfy;
433  tau33(i,j,k) = ( w(i, j, k+1) - w(i, j, k) )*idz0;
434  });
435 
436  Box planexy = tbxxy; planexy.setBig(2, planexy.smallEnd(2) );
437  tbxxy.growLo(2,-1);
438 
439  ParallelFor(planexy,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
440  Real mfy = 0.5 * (mf_uy(i,j,0) + mf_uy(i ,j-1,0));
441  Real mfx = 0.5 * (mf_vx(i,j,0) + mf_vx(i-1,j ,0));
442 
443  tau12(i,j,k) = 0.5 * ( (u(i, j, k) - u(i , j-1, k))*dxInv[1]*mfy
444  + (v(i, j, k) - v(i-1, j , k))*dxInv[0]*mfx);
445  tau21(i,j,k) = tau12(i,j,k);
446  });
447  }
448 
449  //***********************************************************************************
450  // Z-lo w/out Z-Dirichlet
451  //***********************************************************************************
452  if (!zl_u_dir && (tbxxz.smallEnd(2) == domain_xz.smallEnd(2)) ) {
453  Box planexz = tbxxz; planexz.setBig(2, planexz.smallEnd(2));
454  tbxxz.growLo(2,-1);
455 
456  ParallelFor(planexz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
457  Real mfx = mf_ux(i,j,0);
458 
459  Real dz_inv = (k == 0) ? 1.0 / dz_ptr[k] : 2.0 / (dz_ptr[k] + dz_ptr[k-1]);
460 
461  Real du_dz = (u(i, j, k) - u(i , j, k-1))*dz_inv;
462  tau13(i,j,k) = 0.5 * ( du_dz
463  + (w(i, j, k) - w(i-1, j, k ))*dxInv[0] * mfx );
464  tau31(i,j,k) = tau13(i,j,k);
465 
466  if (tau13i) tau13i(i,j,k) = 0.5 * du_dz;
467  });
468  }
469  if (!zl_v_dir && (tbxyz.smallEnd(2) == domain_yz.smallEnd(2))) {
470  Box planeyz = tbxyz; planeyz.setBig(2, planeyz.smallEnd(2) );
471  tbxyz.growLo(2,-1);
472 
473  ParallelFor(planeyz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
474  Real mfy = mf_vy(i,j,0);
475 
476  Real dz_inv = (k == 0) ? 1.0 / dz_ptr[k] : 2.0 / (dz_ptr[k] + dz_ptr[k-1]);
477 
478  Real dv_dz = (v(i, j, k) - v(i, j , k-1))*dz_inv;
479  tau23(i,j,k) = 0.5 * ( dv_dz
480  + (w(i, j, k) - w(i, j-1, k ))*dxInv[1] * mfy );
481  tau32(i,j,k) = tau23(i,j,k);
482 
483  if (tau23i) tau23i(i,j,k) = 0.5 * dv_dz;
484  });
485  }
486 
487  //***********************************************************************************
488  // Z-hi w/out Z-Dirichlet (h_xi = h_eta = 0)
489  //***********************************************************************************
490  if (!zh_u_dir && (tbxxz.bigEnd(2) == domain_xz.bigEnd(2))) {
491  Box planexz = tbxxz; planexz.setSmall(2, planexz.bigEnd(2) );
492  tbxxz.growHi(2,-1);
493 
494  ParallelFor(planexz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
495  Real mfx = mf_ux(i,j,0);
496 
497  Real dz_inv = (k == 0) ? 1.0 / dz_ptr[k] : 2.0 / (dz_ptr[k] + dz_ptr[k-1]);
498 
499  Real du_dz = (u(i, j, k) - u(i , j, k-1))*dz_inv;
500  tau13(i,j,k) = 0.5 * ( du_dz
501  + (w(i, j, k) - w(i-1, j, k ))*dxInv[0]*mfx );
502  tau31(i,j,k) = tau13(i,j,k);
503 
504  if (tau13i) tau13i(i,j,k) = 0.5 * du_dz;
505  });
506  }
507  if (!zh_v_dir && (tbxyz.bigEnd(2) == domain_yz.bigEnd(2))) {
508  Box planeyz = tbxyz; planeyz.setSmall(2, planeyz.bigEnd(2) );
509  tbxyz.growHi(2,-1);
510 
511  ParallelFor(planeyz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
512  Real mfy = mf_vy(i,j,0);
513 
514  Real dz_inv = (k == 0) ? 1.0 / dz_ptr[k] : 2.0 / (dz_ptr[k] + dz_ptr[k-1]);
515 
516  Real dv_dz = (v(i, j, k) - v(i, j , k-1))*dz_inv;
517  tau23(i,j,k) = 0.5 * ( dv_dz
518  + (w(i, j, k) - w(i, j-1, k ))*dxInv[1]*mfy );
519  tau32(i,j,k) = tau23(i,j,k);
520 
521  if (tau23i) tau23i(i,j,k) = 0.5 * dv_dz;
522  });
523  }
524 
525  //***********************************************************************************
526  // Fill the interior cells
527  //***********************************************************************************
528  // Cell centered strains
529  ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
530  Real mfx = mf_mx(i,j,0);
531  Real mfy = mf_my(i,j,0);
532 
533  Real dz_inv = 1.0 / dz_ptr[k];
534 
535  tau11(i,j,k) = (u(i+1, j, k) - u(i, j, k))*dxInv[0] * mfx;
536  tau22(i,j,k) = (v(i, j+1, k) - v(i, j, k))*dxInv[1] * mfy;
537  tau33(i,j,k) = (w(i, j, k+1) - w(i, j, k))*dz_inv;
538  });
539 
540  // Off-diagonal strains
541  ParallelFor(tbxxy,tbxxz,tbxyz,
542  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
543  Real mfy = 0.5 * (mf_uy(i,j,0) + mf_uy(i ,j-1,0));
544  Real mfx = 0.5 * (mf_vx(i,j,0) + mf_vx(i-1,j ,0));
545 
546  tau12(i,j,k) = 0.5 * ( (u(i, j, k) - u(i , j-1, k))*dxInv[1]*mfy
547  + (v(i, j, k) - v(i-1, j , k))*dxInv[0]*mfx);
548  tau21(i,j,k) = tau12(i,j,k);
549  },
550  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
551  Real mfx = mf_ux(i,j,0);
552 
553  Real dz_inv = (k == 0) ? 1.0 / dz_ptr[k] : 2.0 / (dz_ptr[k] + dz_ptr[k-1]);
554 
555  Real du_dz = (u(i, j, k) - u(i , j, k-1))*dz_inv;
556  tau13(i,j,k) = 0.5 * ( du_dz
557  + (w(i, j, k) - w(i-1, j, k ))*dxInv[0] * mfx );
558  tau31(i,j,k) = tau13(i,j,k);
559 
560  if (tau13i) tau13i(i,j,k) = 0.5 * du_dz;
561  },
562  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
563  Real mfy = mf_vy(i,j,0);
564 
565  Real dz_inv = (k == 0) ? 1.0 / dz_ptr[k] : 2.0 / (dz_ptr[k] + dz_ptr[k-1]);
566 
567  Real dv_dz = (v(i, j, k) - v(i, j , k-1))*dz_inv;
568  tau23(i,j,k) = 0.5 * ( dv_dz
569  + (w(i, j, k) - w(i, j-1, k ))*dxInv[1] * mfy );
570  tau32(i,j,k) = tau23(i,j,k);
571 
572  if (tau23i) tau23i(i,j,k) = 0.5 * dv_dz;
573  });
574 }
tau12
@ tau12
Definition: ERF_DataStruct.H:31
tau23
@ tau23
Definition: ERF_DataStruct.H:31
tau33
@ tau33
Definition: ERF_DataStruct.H:31
tau22
@ tau22
Definition: ERF_DataStruct.H:31
tau11
@ tau11
Definition: ERF_DataStruct.H:31
tau32
@ tau32
Definition: ERF_DataStruct.H:31
tau31
@ tau31
Definition: ERF_DataStruct.H:31
tau21
@ tau21
Definition: ERF_DataStruct.H:31
tau13
@ tau13
Definition: ERF_DataStruct.H:31
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 caller graph for this function: