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

Functions

void ComputeStressConsVisc_T (Box bxcc, Box tbxxy, Box tbxxz, Box tbxyz, Real mu_eff, const Array4< const Real > &cell_data, 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 Array4< const Real > &er_arr, const Array4< const Real > &z_nd, const Array4< const Real > &detJ, 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, Array4< Real > &tau13i, Array4< Real > &tau23i, Array4< Real > &tau33i)
 
void ComputeStressVarVisc_T (Box bxcc, Box tbxxy, Box tbxxz, Box tbxyz, Real mu_eff, const Array4< const Real > &mu_turb, const Array4< const Real > &cell_data, 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 Array4< const Real > &er_arr, const Array4< const Real > &z_nd, const Array4< const Real > &detJ, 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, Array4< Real > &tau13i, Array4< Real > &tau23i, Array4< Real > &tau33i)
 

Function Documentation

◆ ComputeStressConsVisc_T()

void ComputeStressConsVisc_T ( Box  bxcc,
Box  tbxxy,
Box  tbxxz,
Box  tbxyz,
Real  mu_eff,
const Array4< const Real > &  cell_data,
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 Array4< const Real > &  er_arr,
const Array4< const Real > &  z_nd,
const Array4< const Real > &  detJ,
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,
Array4< Real > &  tau13i,
Array4< Real > &  tau23i,
Array4< Real > &  tau33i 
)

Function for computing the stress with constant viscosity and 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]mu_effconstant molecular viscosity
[in]cell_datato access rho if ConstantAlpha
[in,out]tau1111 strain -> stress
[in,out]tau2222 strain -> stress
[in,out]tau3333 strain -> stress
[in,out]tau1212 strain -> stress
[in,out]tau1313 strain -> stress
[in,out]tau2121 strain -> stress
[in,out]tau2323 strain -> stress
[in,out]tau3131 strain -> stress
[in,out]tau3232 strain -> stress
[in]er_arrexpansion rate
[in]z_ndnodal array of physical z heights
[in]dxInvinverse cell size array
[in,out]tau13icontribution to stress from du/dz
[in,out]tau23icontribution to stress from dv/dz
[in,out]tau33icontribution to stress from dw/dz
51 {
52  // NOTE: mu_eff includes factor of 2
53 
54  // Handle constant alpha case, in which the provided mu_eff is actually
55  // "alpha" and the viscosity needs to be scaled by rho. This can be further
56  // optimized with if statements below instead of creating a new FAB,
57  // but this is implementation is cleaner.
58  FArrayBox temp;
59  Box gbx = bxcc; // Note: bxcc have been grown in x/y only.
60  gbx.grow(IntVect(0,0,1));
61  temp.resize(gbx,1, The_Async_Arena());
62  Array4<Real> rhoAlpha = temp.array();
63 
64  if (cell_data)
65  // constant alpha (stored in mu_eff)
66  {
67  ParallelFor(gbx,
68  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
69  rhoAlpha(i,j,k) = cell_data(i, j, k, Rho_comp) * mu_eff;
70  });
71  }
72  else
73  // constant mu_eff
74  {
75  ParallelFor(gbx,
76  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
77  rhoAlpha(i,j,k) = mu_eff;
78  });
79  }
80 
81  //***********************************************************************************
82  // NOTE: The first block computes (S-D).
83  // The second block computes 2mu*JT*(S-D)
84  // Boxes are copied here for extrapolations in the second block operations
85  //***********************************************************************************
86  Box bxcc2 = bxcc;
87  bxcc2.grow(IntVect(-1,-1,0));
88 
89  // First block: compute S-D
90  //***********************************************************************************
91  Real OneThird = (1./3.);
92  ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
93  if (tau33i) tau33i(i,j,k) = tau33(i,j,k);
94 
95  tau11(i,j,k) -= OneThird*er_arr(i,j,k);
96  tau22(i,j,k) -= OneThird*er_arr(i,j,k);
97  tau33(i,j,k) -= OneThird*er_arr(i,j,k);
98  });
99 
100  // Second block: compute 2mu*JT*(S-D)
101  //***********************************************************************************
102  // Fill tau33 first (no linear combination extrapolation)
103  //-----------------------------------------------------------------------------------
104  ParallelFor(bxcc2,
105  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
106  {
107  Real mfx = mf_mx(i,j,0);
108  Real mfy = mf_my(i,j,0);
109 
110  Real met_h_xi,met_h_eta;
111  met_h_xi = Compute_h_xi_AtCellCenter (i,j,k,dxInv,z_nd);
112  met_h_eta = Compute_h_eta_AtCellCenter (i,j,k,dxInv,z_nd);
113 
114  Real tau31bar = 0.25 * ( tau31(i , j , k ) + tau31(i+1, j , k )
115  + tau31(i , j , k+1) + tau31(i+1, j , k+1) );
116  Real tau32bar = 0.25 * ( tau32(i , j , k ) + tau32(i , j+1, k )
117  + tau32(i , j , k+1) + tau32(i , j+1, k+1) );
118  Real mu_tot = rhoAlpha(i,j,k);
119 
120  tau33(i,j,k) -= met_h_xi*mfx*tau31bar + met_h_eta*mfy*tau32bar;
121  tau33(i,j,k) *= -mu_tot;
122 
123  if (tau33i) tau33i(i,j,k) *= -mu_tot;
124  });
125 
126  // Second block: compute 2mu*JT*(S-D)
127  //***********************************************************************************
128  // Fill tau13, tau23 next (linear combination extrapolation)
129  //-----------------------------------------------------------------------------------
130  // Extrapolate tau13 & tau23 to bottom
131  {
132  Box planexz = tbxxz; planexz.setBig(2, planexz.smallEnd(2) );
133  tbxxz.growLo(2,-1);
134  ParallelFor(planexz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
135  {
136  Real mfx = mf_ux(i,j,0);
137  Real mfy = mf_uy(i,j,0);
138 
139  Real met_h_xi,met_h_eta,met_h_zeta;
140  met_h_xi = Compute_h_xi_AtEdgeCenterJ (i,j,k,dxInv,z_nd);
141  met_h_eta = Compute_h_eta_AtEdgeCenterJ (i,j,k,dxInv,z_nd);
142  met_h_zeta = Compute_h_zeta_AtEdgeCenterJ(i,j,k,dxInv,z_nd);
143 
144  Real tau11lo = 0.5 * ( tau11(i , j , k ) + tau11(i-1, j , k ) );
145  Real tau11hi = 0.5 * ( tau11(i , j , k+1) + tau11(i-1, j , k+1) );
146  Real tau11bar = 1.5*tau11lo - 0.5*tau11hi;
147 
148  Real tau12lo = 0.5 * ( tau12(i , j , k ) + tau12(i , j+1, k ) );
149  Real tau12hi = 0.5 * ( tau12(i , j , k+1) + tau12(i , j+1, k+1) );
150  Real tau12bar = 1.5*tau12lo - 0.5*tau12hi;
151 
152  Real mu_tot = 0.25*( rhoAlpha(i-1, j, k ) + rhoAlpha(i, j, k )
153  + rhoAlpha(i-1, j, k-1) + rhoAlpha(i, j, k-1) );
154 
155  tau13(i,j,k) -= met_h_xi*mfx*tau11bar + met_h_eta*mfy*tau12bar;
156  tau13(i,j,k) *= -mu_tot;
157  if (tau13i) tau13i(i,j,k) *= -mu_tot;
158 
159  tau31(i,j,k) *= -mu_tot*met_h_zeta/mfy;
160  });
161 
162  Box planeyz = tbxyz; planeyz.setBig(2, planeyz.smallEnd(2) );
163  tbxyz.growLo(2,-1);
164  ParallelFor(planeyz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
165  {
166  Real mfx = mf_vx(i,j,0);
167  Real mfy = mf_vy(i,j,0);
168 
169  Real met_h_xi,met_h_eta,met_h_zeta;
170  met_h_xi = Compute_h_xi_AtEdgeCenterI (i,j,k,dxInv,z_nd);
171  met_h_eta = Compute_h_eta_AtEdgeCenterI (i,j,k,dxInv,z_nd);
172  met_h_zeta = Compute_h_zeta_AtEdgeCenterI(i,j,k,dxInv,z_nd);
173 
174  Real tau21lo = 0.5 * ( tau21(i , j , k ) + tau21(i+1, j , k ) );
175  Real tau21hi = 0.5 * ( tau21(i , j , k+1) + tau21(i+1, j , k+1) );
176  Real tau21bar = 1.5*tau21lo - 0.5*tau21hi;
177 
178  Real tau22lo = 0.5 * ( tau22(i , j , k ) + tau22(i , j-1, k ) );
179  Real tau22hi = 0.5 * ( tau22(i , j , k+1) + tau22(i , j-1, k+1) );
180  Real tau22bar = 1.5*tau22lo - 0.5*tau22hi;
181 
182  Real mu_tot = 0.25*( rhoAlpha(i, j-1, k ) + rhoAlpha(i, j, k )
183  + rhoAlpha(i, j-1, k-1) + rhoAlpha(i, j, k-1) );
184 
185  tau23(i,j,k) -= met_h_xi*mfx*tau21bar + met_h_eta*mfy*tau22bar;
186  tau23(i,j,k) *= -mu_tot;
187  if (tau23i) tau23i(i,j,k) *= -mu_tot;
188 
189  tau32(i,j,k) *= -mu_tot*met_h_zeta/mfx;
190  });
191  }
192  // Extrapolate tau13 & tau23 to top
193  {
194  Box planexz = tbxxz; planexz.setSmall(2, planexz.bigEnd(2) );
195  tbxxz.growHi(2,-1);
196  ParallelFor(planexz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
197  {
198  Real mfx = mf_ux(i,j,0);
199  Real mfy = mf_uy(i,j,0);
200 
201  Real met_h_xi,met_h_eta,met_h_zeta;
202  met_h_xi = Compute_h_xi_AtEdgeCenterJ (i,j,k,dxInv,z_nd);
203  met_h_eta = Compute_h_eta_AtEdgeCenterJ (i,j,k,dxInv,z_nd);
204  met_h_zeta = Compute_h_zeta_AtEdgeCenterJ(i,j,k,dxInv,z_nd);
205 
206  Real tau11lo = 0.5 * ( tau11(i , j , k-2) + tau11(i-1, j , k-2) );
207  Real tau11hi = 0.5 * ( tau11(i , j , k-1) + tau11(i-1, j , k-1) );
208  Real tau11bar = 1.5*tau11hi - 0.5*tau11lo;
209 
210  Real tau12lo = 0.5 * ( tau12(i , j , k-2) + tau12(i , j+1, k-2) );
211  Real tau12hi = 0.5 * ( tau12(i , j , k-1) + tau12(i , j+1, k-1) );
212  Real tau12bar = 1.5*tau12hi - 0.5*tau12lo;
213 
214  Real mu_tot = 0.25*( rhoAlpha(i-1, j, k ) + rhoAlpha(i, j, k )
215  + rhoAlpha(i-1, j, k-1) + rhoAlpha(i, j, k-1) );
216 
217  tau13(i,j,k) -= met_h_xi*mfx*tau11bar + met_h_eta*mfy*tau12bar;
218  tau13(i,j,k) *= -mu_tot;
219  if (tau13i) tau13i(i,j,k) *= -mu_tot;
220 
221  tau31(i,j,k) *= -mu_tot*met_h_zeta/mfy;
222  });
223 
224  Box planeyz = tbxyz; planeyz.setSmall(2, planeyz.bigEnd(2) );
225  tbxyz.growHi(2,-1);
226  ParallelFor(planeyz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
227  {
228  Real mfx = mf_vx(i,j,0);
229  Real mfy = mf_vy(i,j,0);
230 
231  Real met_h_xi,met_h_eta,met_h_zeta;
232  met_h_xi = Compute_h_xi_AtEdgeCenterI (i,j,k,dxInv,z_nd);
233  met_h_eta = Compute_h_eta_AtEdgeCenterI (i,j,k,dxInv,z_nd);
234  met_h_zeta = Compute_h_zeta_AtEdgeCenterI(i,j,k,dxInv,z_nd);
235 
236  Real tau21lo = 0.5 * ( tau21(i , j , k-2) + tau21(i+1, j , k-2) );
237  Real tau21hi = 0.5 * ( tau21(i , j , k-1) + tau21(i+1, j , k-1) );
238  Real tau21bar = 1.5*tau21hi - 0.5*tau21lo;
239 
240  Real tau22lo = 0.5 * ( tau22(i , j , k-2) + tau22(i , j-1, k-2) );
241  Real tau22hi = 0.5 * ( tau22(i , j , k-1) + tau22(i , j-1, k-1) );
242  Real tau22bar = 1.5*tau22hi - 0.5*tau22lo;
243 
244  Real mu_tot = 0.25*( rhoAlpha(i, j-1, k ) + rhoAlpha(i, j, k )
245  + rhoAlpha(i, j-1, k-1) + rhoAlpha(i, j, k-1) );
246 
247  tau23(i,j,k) -= met_h_xi*mfx*tau21bar + met_h_eta*mfy*tau22bar;
248  tau23(i,j,k) *= -mu_tot;
249  if (tau23i) tau23i(i,j,k) *= -mu_tot;
250 
251  tau32(i,j,k) *= -mu_tot*met_h_zeta/mfx;
252  });
253  }
254 
255  // Second block: compute 2mu*JT*(S-D)
256  //***********************************************************************************
257  // Fill tau13, tau23 next (valid averaging region)
258  //-----------------------------------------------------------------------------------
259  ParallelFor(tbxxz,tbxyz,
260  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
261  {
262  Real mfx = mf_ux(i,j,0);
263  Real mfy = mf_uy(i,j,0);
264 
265  Real met_h_xi,met_h_eta,met_h_zeta;
266  met_h_xi = Compute_h_xi_AtEdgeCenterJ (i,j,k,dxInv,z_nd);
267  met_h_eta = Compute_h_eta_AtEdgeCenterJ (i,j,k,dxInv,z_nd);
268  met_h_zeta = Compute_h_zeta_AtEdgeCenterJ(i,j,k,dxInv,z_nd);
269 
270  Real tau11bar = 0.25 * ( tau11(i , j , k ) + tau11(i-1, j , k )
271  + tau11(i , j , k-1) + tau11(i-1, j , k-1) );
272  Real tau12bar = 0.25 * ( tau12(i , j , k ) + tau12(i , j+1, k )
273  + tau12(i , j , k-1) + tau12(i , j+1, k-1) );
274  Real mu_tot = 0.25 * ( rhoAlpha(i-1, j , k ) + rhoAlpha(i , j , k )
275  + rhoAlpha(i-1, j , k-1) + rhoAlpha(i , j , k-1) );
276 
277  tau13(i,j,k) -= met_h_xi*mfx*tau11bar + met_h_eta*mfy*tau12bar;
278  tau13(i,j,k) *= -mu_tot;
279  if (tau13i) tau13i(i,j,k) *= -mu_tot;
280 
281  tau31(i,j,k) *= -mu_tot*met_h_zeta/mfy;
282  },
283  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
284  {
285  Real mfx = mf_vx(i,j,0);
286  Real mfy = mf_vy(i,j,0);
287 
288  Real met_h_xi,met_h_eta,met_h_zeta;
289  met_h_xi = Compute_h_xi_AtEdgeCenterI (i,j,k,dxInv,z_nd);
290  met_h_eta = Compute_h_eta_AtEdgeCenterI (i,j,k,dxInv,z_nd);
291  met_h_zeta = Compute_h_zeta_AtEdgeCenterI(i,j,k,dxInv,z_nd);
292 
293  Real tau21bar = 0.25 * ( tau21(i , j , k ) + tau21(i+1, j , k )
294  + tau21(i , j , k-1) + tau21(i+1, j , k-1) );
295  Real tau22bar = 0.25 * ( tau22(i , j , k ) + tau22(i , j-1, k )
296  + tau22(i , j , k-1) + tau22(i , j-1, k-1) );
297  Real mu_tot = 0.25 * ( rhoAlpha(i , j-1, k ) + rhoAlpha(i , j , k )
298  + rhoAlpha(i , j-1, k-1) + rhoAlpha(i , j , k-1) );
299 
300  tau23(i,j,k) -= met_h_xi*mfx*tau21bar + met_h_eta*mfy*tau22bar;
301  tau23(i,j,k) *= -mu_tot;
302  if (tau23i) tau23i(i,j,k) *= -mu_tot;
303 
304  tau32(i,j,k) *= -mu_tot*met_h_zeta/mfx;
305  });
306 
307  // Fill the remaining components: tau11, tau22, tau12/21
308  //-----------------------------------------------------------------------------------
309  ParallelFor(bxcc,tbxxy,
310  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
311  {
312  Real mfx = mf_mx(i,j,0);
313  Real mfy = mf_my(i,j,0);
314 
315  Real met_h_zeta = detJ(i,j,k);
316  Real mu_tot = rhoAlpha(i,j,k);
317 
318  tau11(i,j,k) *= -mu_tot*met_h_zeta/mfy;
319  tau22(i,j,k) *= -mu_tot*met_h_zeta/mfx;
320  },
321  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
322  {
323  Real mfx = 0.5 * (mf_ux(i,j,0) + mf_ux(i,j-1,0));
324  Real mfy = 0.5 * (mf_vy(i,j,0) + mf_vy(i-1,j,0));
325 
326  Real met_h_zeta = Compute_h_zeta_AtEdgeCenterK(i,j,k,dxInv,z_nd);
327 
328  Real mu_tot = 0.25*( rhoAlpha(i-1, j , k) + rhoAlpha(i, j , k)
329  + rhoAlpha(i-1, j-1, k) + rhoAlpha(i, j-1, k) );
330 
331  tau12(i,j,k) *= -mu_tot*met_h_zeta/mfx;
332  tau21(i,j,k) *= -mu_tot*met_h_zeta/mfy;
333  });
334 }
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
Rho_comp
#define Rho_comp
Definition: ERF_IndexDefines.H:36
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:19
Compute_h_eta_AtEdgeCenterJ
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_eta_AtEdgeCenterJ(const int &i, const int &j, const int &k, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &z_nd)
Definition: ERF_TerrainMetrics.H:302
Compute_h_xi_AtEdgeCenterJ
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_xi_AtEdgeCenterJ(const int &i, const int &j, const int &k, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &z_nd)
Definition: ERF_TerrainMetrics.H:287
Compute_h_zeta_AtEdgeCenterJ
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_zeta_AtEdgeCenterJ(const int &i, const int &j, const int &k, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &z_nd)
Definition: ERF_TerrainMetrics.H:273
Compute_h_zeta_AtEdgeCenterK
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_zeta_AtEdgeCenterK(const int &i, const int &j, const int &k, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &z_nd)
Definition: ERF_TerrainMetrics.H:228
Compute_h_zeta_AtEdgeCenterI
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_zeta_AtEdgeCenterI(const int &i, const int &j, const int &k, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &z_nd)
Definition: ERF_TerrainMetrics.H:317
Compute_h_eta_AtCellCenter
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_eta_AtCellCenter(const int &i, const int &j, const int &k, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &z_nd)
Definition: ERF_TerrainMetrics.H:83
Compute_h_xi_AtEdgeCenterI
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_xi_AtEdgeCenterI(const int &i, const int &j, const int &k, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &z_nd)
Definition: ERF_TerrainMetrics.H:331
Compute_h_xi_AtCellCenter
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_xi_AtCellCenter(const int &i, const int &j, const int &k, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &z_nd)
Definition: ERF_TerrainMetrics.H:68
Compute_h_eta_AtEdgeCenterI
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_eta_AtEdgeCenterI(const int &i, const int &j, const int &k, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &z_nd)
Definition: ERF_TerrainMetrics.H:344

Referenced by erf_make_tau_terms().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ ComputeStressVarVisc_T()

void ComputeStressVarVisc_T ( Box  bxcc,
Box  tbxxy,
Box  tbxxz,
Box  tbxyz,
Real  mu_eff,
const Array4< const Real > &  mu_turb,
const Array4< const Real > &  cell_data,
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 Array4< const Real > &  er_arr,
const Array4< const Real > &  z_nd,
const Array4< const Real > &  detJ,
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,
Array4< Real > &  tau13i,
Array4< Real > &  tau23i,
Array4< Real > &  tau33i 
)

Function for computing the stress with constant viscosity and 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]mu_effconstant molecular viscosity
[in]mu_turbvariable turbulent viscosity
[in]cell_datato access rho if ConstantAlpha
[in,out]tau1111 strain -> stress
[in,out]tau2222 strain -> stress
[in,out]tau3333 strain -> stress
[in,out]tau1212 strain -> stress
[in,out]tau1313 strain -> stress
[in,out]tau2121 strain -> stress
[in,out]tau2323 strain -> stress
[in,out]tau3131 strain -> stress
[in,out]tau3232 strain -> stress
[in]er_arrexpansion rate
[in]z_ndnodal array of physical z heights
[in]dxInvinverse cell size array
[in,out]tau13icontribution to stress from du/dz
[in,out]tau23icontribution to stress from dv/dz
[in,out]tau33icontribution to stress from dw/dz
383 {
384  // NOTE: mu_eff includes factor of 2
385 
386  // Handle constant alpha case, in which the provided mu_eff is actually
387  // "alpha" and the viscosity needs to be scaled by rho. This can be further
388  // optimized with if statements below instead of creating a new FAB,
389  // but this is implementation is cleaner.
390  FArrayBox temp;
391  Box gbx = bxcc; // Note: bxcc have been grown in x/y only.
392  gbx.grow(IntVect(0,0,1));
393  temp.resize(gbx,1, The_Async_Arena());
394  Array4<Real> rhoAlpha = temp.array();
395 
396  if (cell_data)
397  // constant alpha (stored in mu_eff)
398  {
399  ParallelFor(gbx,
400  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
401  rhoAlpha(i,j,k) = cell_data(i, j, k, Rho_comp) * mu_eff;
402  });
403  }
404  else
405  // constant mu_eff
406  {
407  ParallelFor(gbx,
408  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
409  rhoAlpha(i,j,k) = mu_eff;
410  });
411  }
412 
413  //***********************************************************************************
414  // NOTE: The first block computes (S-D).
415  // The second block computes 2mu*JT*(S-D)
416  // Boxes are copied here for extrapolations in the second block operations
417  //***********************************************************************************
418  Box bxcc2 = bxcc;
419  bxcc2.grow(IntVect(-1,-1,0));
420 
421  // First block: compute S-D
422  //***********************************************************************************
423  Real OneThird = (1./3.);
424  ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
425  if (tau33i) tau33i(i,j,k) = tau33(i,j,k);
426 
427  tau11(i,j,k) -= OneThird*er_arr(i,j,k);
428  tau22(i,j,k) -= OneThird*er_arr(i,j,k);
429  tau33(i,j,k) -= OneThird*er_arr(i,j,k);
430  });
431 
432  // Second block: compute 2mu*JT*(S-D)
433  //***********************************************************************************
434  // Fill tau33 first (no linear combination extrapolation)
435  //-----------------------------------------------------------------------------------
436  ParallelFor(bxcc2,
437  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
438  {
439  Real mfx = mf_mx(i,j,0);
440  Real mfy = mf_my(i,j,0);
441 
442  Real met_h_xi,met_h_eta;
443  met_h_xi = Compute_h_xi_AtCellCenter (i,j,k,dxInv,z_nd);
444  met_h_eta = Compute_h_eta_AtCellCenter (i,j,k,dxInv,z_nd);
445 
446  Real tau31bar = 0.25 * ( tau31(i , j , k ) + tau31(i+1, j , k )
447  + tau31(i , j , k+1) + tau31(i+1, j , k+1) );
448  Real tau32bar = 0.25 * ( tau32(i , j , k ) + tau32(i , j+1, k )
449  + tau32(i , j , k+1) + tau32(i , j+1, k+1) );
450 
451  Real mu_tot = rhoAlpha(i,j,k) + 2.0*mu_turb(i, j, k, EddyDiff::Mom_v);
452 
453  tau33(i,j,k) -= met_h_xi*mfx*tau31bar + met_h_eta*mfy*tau32bar;
454  tau33(i,j,k) *= -mu_tot;
455 
456  if (tau33i) tau33i(i,j,k) *= -mu_tot;
457  });
458 
459  // Second block: compute 2mu*JT*(S-D)
460  //***********************************************************************************
461  // Fill tau13, tau23 next (linear combination extrapolation)
462  //-----------------------------------------------------------------------------------
463  // Extrapolate tau13 & tau23 to bottom
464  {
465  Box planexz = tbxxz; planexz.setBig(2, planexz.smallEnd(2) );
466  tbxxz.growLo(2,-1);
467  ParallelFor(planexz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
468  {
469  Real mfx = mf_ux(i,j,0);
470  Real mfy = mf_uy(i,j,0);
471 
472  Real met_h_xi,met_h_eta,met_h_zeta;
473  met_h_xi = Compute_h_xi_AtEdgeCenterJ (i,j,k,dxInv,z_nd);
474  met_h_eta = Compute_h_eta_AtEdgeCenterJ (i,j,k,dxInv,z_nd);
475  met_h_zeta = Compute_h_zeta_AtEdgeCenterJ(i,j,k,dxInv,z_nd);
476 
477  Real tau11lo = 0.5 * ( tau11(i , j , k ) + tau11(i-1, j , k ) );
478  Real tau11hi = 0.5 * ( tau11(i , j , k+1) + tau11(i-1, j , k+1) );
479  Real tau11bar = 1.5*tau11lo - 0.5*tau11hi;
480 
481  Real tau12lo = 0.5 * ( tau12(i , j , k ) + tau12(i , j+1, k ) );
482  Real tau12hi = 0.5 * ( tau12(i , j , k+1) + tau12(i , j+1, k+1) );
483  Real tau12bar = 1.5*tau12lo - 0.5*tau12hi;
484 
485  Real mu_bar = 0.25*( mu_turb(i-1, j, k , EddyDiff::Mom_v) + mu_turb(i, j, k , EddyDiff::Mom_v)
486  + mu_turb(i-1, j, k-1, EddyDiff::Mom_v) + mu_turb(i, j, k-1, EddyDiff::Mom_v) );
487  Real rhoAlpha_bar = 0.25*( rhoAlpha(i-1, j, k ) + rhoAlpha(i, j, k )
488  + rhoAlpha(i-1, j, k-1) + rhoAlpha(i, j, k-1) );
489  Real mu_tot = rhoAlpha_bar + 2.0*mu_bar;
490 
491  tau13(i,j,k) -= met_h_xi*mfx*tau11bar + met_h_eta*mfy*tau12bar;
492  tau13(i,j,k) *= -mu_tot;
493  if (tau13i) tau13i(i,j,k) *= -mu_tot;
494 
495  tau31(i,j,k) *= -mu_tot*met_h_zeta/mfy;
496  });
497 
498  Box planeyz = tbxyz; planeyz.setBig(2, planeyz.smallEnd(2) );
499  tbxyz.growLo(2,-1);
500  ParallelFor(planeyz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
501  {
502  Real mfx = mf_vx(i,j,0);
503  Real mfy = mf_vy(i,j,0);
504 
505  Real met_h_xi,met_h_eta,met_h_zeta;
506  met_h_xi = Compute_h_xi_AtEdgeCenterI (i,j,k,dxInv,z_nd);
507  met_h_eta = Compute_h_eta_AtEdgeCenterI (i,j,k,dxInv,z_nd);
508  met_h_zeta = Compute_h_zeta_AtEdgeCenterI(i,j,k,dxInv,z_nd);
509 
510  Real tau21lo = 0.5 * ( tau21(i , j , k ) + tau21(i+1, j , k ) );
511  Real tau21hi = 0.5 * ( tau21(i , j , k+1) + tau21(i+1, j , k+1) );
512  Real tau21bar = 1.5*tau21lo - 0.5*tau21hi;
513 
514  Real tau22lo = 0.5 * ( tau22(i , j , k ) + tau22(i , j-1, k ) );
515  Real tau22hi = 0.5 * ( tau22(i , j , k+1) + tau22(i , j-1, k+1) );
516  Real tau22bar = 1.5*tau22lo - 0.5*tau22hi;
517 
518  Real mu_bar = 0.25*( mu_turb(i, j-1, k , EddyDiff::Mom_v) + mu_turb(i, j, k , EddyDiff::Mom_v)
519  + mu_turb(i, j-1, k-1, EddyDiff::Mom_v) + mu_turb(i, j, k-1, EddyDiff::Mom_v) );
520  Real rhoAlpha_bar = 0.25*( rhoAlpha(i, j-1, k ) + rhoAlpha(i, j, k )
521  + rhoAlpha(i, j-1, k-1) + rhoAlpha(i, j, k-1) );
522  Real mu_tot = rhoAlpha_bar + 2.0*mu_bar;
523 
524  tau23(i,j,k) -= met_h_xi*mfx*tau21bar + met_h_eta*mfy*tau22bar;
525  tau23(i,j,k) *= -mu_tot;
526  if (tau23i) tau23i(i,j,k) *= -mu_tot;
527 
528  tau32(i,j,k) *= -mu_tot*met_h_zeta/mfx;
529  });
530  }
531  // Extrapolate tau13 & tau23 to top
532  {
533  Box planexz = tbxxz; planexz.setSmall(2, planexz.bigEnd(2) );
534  tbxxz.growHi(2,-1);
535  ParallelFor(planexz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
536  {
537  Real mfx = mf_ux(i,j,0);
538  Real mfy = mf_uy(i,j,0);
539 
540  Real met_h_xi,met_h_eta,met_h_zeta;
541  met_h_xi = Compute_h_xi_AtEdgeCenterJ (i,j,k,dxInv,z_nd);
542  met_h_eta = Compute_h_eta_AtEdgeCenterJ (i,j,k,dxInv,z_nd);
543  met_h_zeta = Compute_h_zeta_AtEdgeCenterJ(i,j,k,dxInv,z_nd);
544 
545  Real tau11lo = 0.5 * ( tau11(i , j , k-2) + tau11(i-1, j , k-2) );
546  Real tau11hi = 0.5 * ( tau11(i , j , k-1) + tau11(i-1, j , k-1) );
547  Real tau11bar = 1.5*tau11hi - 0.5*tau11lo;
548 
549  Real tau12lo = 0.5 * ( tau12(i , j , k-2) + tau12(i , j+1, k-2) );
550  Real tau12hi = 0.5 * ( tau12(i , j , k-1) + tau12(i , j+1, k-1) );
551  Real tau12bar = 1.5*tau12hi - 0.5*tau12lo;
552 
553  Real mu_bar = 0.25*( mu_turb(i-1, j, k , EddyDiff::Mom_v) + mu_turb(i, j, k , EddyDiff::Mom_v)
554  + mu_turb(i-1, j, k-1, EddyDiff::Mom_v) + mu_turb(i, j, k-1, EddyDiff::Mom_v) );
555  Real rhoAlpha_bar = 0.25*( rhoAlpha(i-1, j, k ) + rhoAlpha(i, j, k )
556  + rhoAlpha(i-1, j, k-1) + rhoAlpha(i, j, k-1) );
557  Real mu_tot = rhoAlpha_bar + 2.0*mu_bar;
558 
559  tau13(i,j,k) -= met_h_xi*mfx*tau11bar + met_h_eta*mfy*tau12bar;
560  tau13(i,j,k) *= -mu_tot;
561  if (tau13i) tau13i(i,j,k) *= -mu_tot;
562 
563  tau31(i,j,k) *= -mu_tot*met_h_zeta/mfy;
564  });
565 
566  Box planeyz = tbxyz; planeyz.setSmall(2, planeyz.bigEnd(2) );
567  tbxyz.growHi(2,-1);
568  ParallelFor(planeyz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
569  {
570  Real mfx = mf_vx(i,j,0);
571  Real mfy = mf_vy(i,j,0);
572 
573  Real met_h_xi,met_h_eta,met_h_zeta;
574  met_h_xi = Compute_h_xi_AtEdgeCenterI (i,j,k,dxInv,z_nd);
575  met_h_eta = Compute_h_eta_AtEdgeCenterI (i,j,k,dxInv,z_nd);
576  met_h_zeta = Compute_h_zeta_AtEdgeCenterI(i,j,k,dxInv,z_nd);
577 
578  Real tau21lo = 0.5 * ( tau21(i , j , k-2) + tau21(i+1, j , k-2) );
579  Real tau21hi = 0.5 * ( tau21(i , j , k-1) + tau21(i+1, j , k-1) );
580  Real tau21bar = 1.5*tau21hi - 0.5*tau21lo;
581 
582  Real tau22lo = 0.5 * ( tau22(i , j , k-2) + tau22(i , j-1, k-2) );
583  Real tau22hi = 0.5 * ( tau22(i , j , k-1) + tau22(i , j-1, k-1) );
584  Real tau22bar = 1.5*tau22hi - 0.5*tau22lo;
585 
586  Real mu_bar = 0.25*( mu_turb(i, j-1, k , EddyDiff::Mom_v) + mu_turb(i, j, k , EddyDiff::Mom_v)
587  + mu_turb(i, j-1, k-1, EddyDiff::Mom_v) + mu_turb(i, j, k-1, EddyDiff::Mom_v) );
588  Real rhoAlpha_bar = 0.25*( rhoAlpha(i, j-1, k ) + rhoAlpha(i, j, k )
589  + rhoAlpha(i, j-1, k-1) + rhoAlpha(i, j, k-1) );
590  Real mu_tot = rhoAlpha_bar + 2.0*mu_bar;
591 
592  tau23(i,j,k) -= met_h_xi*mfx*tau21bar + met_h_eta*mfy*tau22bar;
593  tau23(i,j,k) *= -mu_tot;
594  if (tau23i) tau23i(i,j,k) *= -mu_tot;
595 
596  tau32(i,j,k) *= -mu_tot*met_h_zeta/mfx;
597  });
598  }
599 
600  // Second block: compute 2mu*JT*(S-D)
601  //***********************************************************************************
602  // Fill tau13, tau23 next (valid averaging region)
603  //-----------------------------------------------------------------------------------
604  ParallelFor(tbxxz,tbxyz,
605  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
606  {
607  Real mfx = mf_ux(i,j,0);
608  Real mfy = mf_uy(i,j,0);
609 
610  Real met_h_xi,met_h_eta,met_h_zeta;
611  met_h_xi = Compute_h_xi_AtEdgeCenterJ (i,j,k,dxInv,z_nd);
612  met_h_eta = Compute_h_eta_AtEdgeCenterJ (i,j,k,dxInv,z_nd);
613  met_h_zeta = Compute_h_zeta_AtEdgeCenterJ(i,j,k,dxInv,z_nd);
614 
615  Real tau11bar = 0.25 * ( tau11(i , j , k ) + tau11(i-1, j , k )
616  + tau11(i , j , k-1) + tau11(i-1, j , k-1) );
617  Real tau12bar = 0.25 * ( tau12(i , j , k ) + tau12(i , j+1, k )
618  + tau12(i , j , k-1) + tau12(i , j+1, k-1) );
619 
620  Real mu_bar = 0.25 * ( mu_turb(i-1, j , k , EddyDiff::Mom_v) + mu_turb(i , j , k , EddyDiff::Mom_v)
621  + mu_turb(i-1, j , k-1, EddyDiff::Mom_v) + mu_turb(i , j , k-1, EddyDiff::Mom_v) );
622  Real rhoAlpha_bar = 0.25 * ( rhoAlpha(i-1, j , k ) + rhoAlpha(i , j , k )
623  + rhoAlpha(i-1, j , k-1) + rhoAlpha(i , j , k-1) );
624  Real mu_tot = rhoAlpha_bar + 2.0*mu_bar;
625 
626  tau13(i,j,k) -= met_h_xi*mfx*tau11bar + met_h_eta*mfy*tau12bar;
627  tau13(i,j,k) *= -mu_tot;
628  if (tau13i) tau13i(i,j,k) *= -mu_tot;
629 
630  tau31(i,j,k) *= -mu_tot*met_h_zeta/mfy;
631  },
632  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
633  {
634  Real mfx = mf_vx(i,j,0);
635  Real mfy = mf_vy(i,j,0);
636 
637  Real met_h_xi,met_h_eta,met_h_zeta;
638  met_h_xi = Compute_h_xi_AtEdgeCenterI (i,j,k,dxInv,z_nd);
639  met_h_eta = Compute_h_eta_AtEdgeCenterI (i,j,k,dxInv,z_nd);
640  met_h_zeta = Compute_h_zeta_AtEdgeCenterI(i,j,k,dxInv,z_nd);
641 
642  Real tau21bar = 0.25 * ( tau21(i , j , k ) + tau21(i+1, j , k )
643  + tau21(i , j , k-1) + tau21(i+1, j , k-1) );
644  Real tau22bar = 0.25 * ( tau22(i , j , k ) + tau22(i , j-1, k )
645  + tau22(i , j , k-1) + tau22(i , j-1, k-1) );
646 
647  Real mu_bar = 0.25 * ( mu_turb(i , j-1, k , EddyDiff::Mom_v) + mu_turb(i , j , k , EddyDiff::Mom_v)
648  + mu_turb(i , j-1, k-1, EddyDiff::Mom_v) + mu_turb(i , j , k-1, EddyDiff::Mom_v) );
649  Real rhoAlpha_bar = 0.25 * ( rhoAlpha(i , j-1, k ) + rhoAlpha(i , j , k )
650  + rhoAlpha(i , j-1, k-1) + rhoAlpha(i , j , k-1) );
651  Real mu_tot = rhoAlpha_bar + 2.0*mu_bar;
652 
653  tau23(i,j,k) -= met_h_xi*mfx*tau21bar + met_h_eta*mfy*tau22bar;
654  tau23(i,j,k) *= -mu_tot;
655  if (tau23i) tau23i(i,j,k) *= -mu_tot;
656 
657  tau32(i,j,k) *= -mu_tot*met_h_zeta/mfx;
658  });
659 
660  // Fill the remaining components: tau11, tau22, tau12/21
661  //-----------------------------------------------------------------------------------
662  ParallelFor(bxcc,tbxxy,
663  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
664  {
665  Real mfx = mf_mx(i,j,0);
666  Real mfy = mf_my(i,j,0);
667 
668  Real met_h_zeta = detJ(i,j,k);
669 
670  Real mu_tot = rhoAlpha(i,j,k) + 2.0*mu_turb(i, j, k, EddyDiff::Mom_h);
671 
672  tau11(i,j,k) *= -mu_tot*met_h_zeta/mfy;
673  tau22(i,j,k) *= -mu_tot*met_h_zeta/mfx;
674  },
675  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
676  {
677  Real mfx = 0.5 * (mf_ux(i,j,0) + mf_ux(i,j-1,0));
678  Real mfy = 0.5 * (mf_vy(i,j,0) + mf_vy(i-1,j,0));
679 
680  Real met_h_zeta = Compute_h_zeta_AtEdgeCenterK(i,j,k,dxInv,z_nd);
681 
682  Real mu_bar = 0.25*( mu_turb(i-1, j , k, EddyDiff::Mom_h) + mu_turb(i, j , k, EddyDiff::Mom_h)
683  + mu_turb(i-1, j-1, k, EddyDiff::Mom_h) + mu_turb(i, j-1, k, EddyDiff::Mom_h) );
684  Real rhoAlpha_bar = 0.25*( rhoAlpha(i-1, j , k) + rhoAlpha(i, j , k)
685  + rhoAlpha(i-1, j-1, k) + rhoAlpha(i, j-1, k) );
686  Real mu_tot = rhoAlpha_bar + 2.0*mu_bar;
687 
688  tau12(i,j,k) *= -mu_tot*met_h_zeta/mfx;
689  tau21(i,j,k) *= -mu_tot*met_h_zeta/mfy;
690  });
691 }
EddyDiff::Mom_h
@ Mom_h
Definition: ERF_IndexDefines.H:170
EddyDiff::Mom_v
@ Mom_v
Definition: ERF_IndexDefines.H:175

Referenced by erf_make_tau_terms().

Here is the call graph for this function:
Here is the caller graph for this function: