ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
Functions
ERF_ComputeStrain_S.cpp File Reference
#include <ERF_Diffusion.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)
 

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 
)

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

Referenced by ERF::advance_dycore(), and erf_make_tau_terms().

Here is the caller graph for this function: