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

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

Referenced by erf_make_tau_terms().

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