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 on a stretched grid.

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