ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
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)
 
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)
 

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 
)

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

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 
)

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
353 {
354  // Handle constant alpha case, in which the provided mu_eff is actually
355  // "alpha" and the viscosity needs to be scaled by rho. This can be further
356  // optimized with if statements below instead of creating a new FAB,
357  // but this is implementation is cleaner.
358  FArrayBox temp;
359  Box gbx = bxcc; // Note: bxcc have been grown in x/y only.
360  gbx.grow(IntVect(0,0,1));
361  temp.resize(gbx,1, The_Async_Arena());
362  Array4<Real> rhoAlpha = temp.array();
363  if (cell_data) {
364  ParallelFor(gbx,
365  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
366  rhoAlpha(i,j,k) = cell_data(i, j, k, Rho_comp) * mu_eff;
367  });
368  } else {
369  ParallelFor(gbx,
370  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
371  rhoAlpha(i,j,k) = mu_eff;
372  });
373  }
374 
375  //***********************************************************************************
376  // NOTE: The first block computes (S-D).
377  // The second block computes 2mu*JT*(S-D)
378  // Boxes are copied here for extrapolations in the second block operations
379  //***********************************************************************************
380  Box bxcc2 = bxcc;
381  bxcc2.grow(IntVect(-1,-1,0));
382 
383  // First block: compute S-D
384  //***********************************************************************************
385  Real OneThird = (1./3.);
386  ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
387  tau11(i,j,k) -= OneThird*er_arr(i,j,k);
388  tau22(i,j,k) -= OneThird*er_arr(i,j,k);
389  tau33(i,j,k) -= OneThird*er_arr(i,j,k);
390  });
391 
392  // Second block: compute 2mu*JT*(S-D)
393  //***********************************************************************************
394  // Fill tau33 first (no linear combination extrapolation)
395  //-----------------------------------------------------------------------------------
396  ParallelFor(bxcc2,
397  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
398  {
399  Real mfx = mf_mx(i,j,0);
400  Real mfy = mf_my(i,j,0);
401 
402  Real met_h_xi,met_h_eta;
403  met_h_xi = Compute_h_xi_AtCellCenter (i,j,k,dxInv,z_nd);
404  met_h_eta = Compute_h_eta_AtCellCenter (i,j,k,dxInv,z_nd);
405 
406  Real tau31bar = 0.25 * ( tau31(i , j , k ) + tau31(i+1, j , k )
407  + tau31(i , j , k+1) + tau31(i+1, j , k+1) );
408  Real tau32bar = 0.25 * ( tau32(i , j , k ) + tau32(i , j+1, k )
409  + tau32(i , j , k+1) + tau32(i , j+1, k+1) );
410 
411  Real mu_tot = rhoAlpha(i,j,k) + 2.0*mu_turb(i, j, k, EddyDiff::Mom_v);
412 
413  tau33(i,j,k) -= met_h_xi*mfx*tau31bar + met_h_eta*mfy*tau32bar;
414  tau33(i,j,k) *= -mu_tot;
415  });
416 
417  // Second block: compute 2mu*JT*(S-D)
418  //***********************************************************************************
419  // Fill tau13, tau23 next (linear combination extrapolation)
420  //-----------------------------------------------------------------------------------
421  // Extrapolate tau13 & tau23 to bottom
422  {
423  Box planexz = tbxxz; planexz.setBig(2, planexz.smallEnd(2) );
424  tbxxz.growLo(2,-1);
425  ParallelFor(planexz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
426  {
427  Real mfx = mf_ux(i,j,0);
428  Real mfy = mf_uy(i,j,0);
429 
430  Real met_h_xi,met_h_eta,met_h_zeta;
431  met_h_xi = Compute_h_xi_AtEdgeCenterJ (i,j,k,dxInv,z_nd);
432  met_h_eta = Compute_h_eta_AtEdgeCenterJ (i,j,k,dxInv,z_nd);
433  met_h_zeta = Compute_h_zeta_AtEdgeCenterJ(i,j,k,dxInv,z_nd);
434 
435  Real tau11lo = 0.5 * ( tau11(i , j , k ) + tau11(i-1, j , k ) );
436  Real tau11hi = 0.5 * ( tau11(i , j , k+1) + tau11(i-1, j , k+1) );
437  Real tau11bar = 1.5*tau11lo - 0.5*tau11hi;
438 
439  Real tau12lo = 0.5 * ( tau12(i , j , k ) + tau12(i , j+1, k ) );
440  Real tau12hi = 0.5 * ( tau12(i , j , k+1) + tau12(i , j+1, k+1) );
441  Real tau12bar = 1.5*tau12lo - 0.5*tau12hi;
442 
443  Real mu_bar = 0.25*( mu_turb(i-1, j, k , EddyDiff::Mom_v) + mu_turb(i, j, k , EddyDiff::Mom_v)
444  + mu_turb(i-1, j, k-1, EddyDiff::Mom_v) + mu_turb(i, j, k-1, EddyDiff::Mom_v) );
445  Real rhoAlpha_bar = 0.25*( rhoAlpha(i-1, j, k ) + rhoAlpha(i, j, k )
446  + rhoAlpha(i-1, j, k-1) + rhoAlpha(i, j, k-1) );
447  Real mu_tot = rhoAlpha_bar + 2.0*mu_bar;
448 
449  tau13(i,j,k) -= met_h_xi*mfx*tau11bar + met_h_eta*mfy*tau12bar;
450  tau13(i,j,k) *= -mu_tot;
451 
452  tau31(i,j,k) *= -mu_tot*met_h_zeta/mfy;
453  });
454 
455  Box planeyz = tbxyz; planeyz.setBig(2, planeyz.smallEnd(2) );
456  tbxyz.growLo(2,-1);
457  ParallelFor(planeyz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
458  {
459  Real mfx = mf_vx(i,j,0);
460  Real mfy = mf_vy(i,j,0);
461 
462  Real met_h_xi,met_h_eta,met_h_zeta;
463  met_h_xi = Compute_h_xi_AtEdgeCenterI (i,j,k,dxInv,z_nd);
464  met_h_eta = Compute_h_eta_AtEdgeCenterI (i,j,k,dxInv,z_nd);
465  met_h_zeta = Compute_h_zeta_AtEdgeCenterI(i,j,k,dxInv,z_nd);
466 
467  Real tau21lo = 0.5 * ( tau21(i , j , k ) + tau21(i+1, j , k ) );
468  Real tau21hi = 0.5 * ( tau21(i , j , k+1) + tau21(i+1, j , k+1) );
469  Real tau21bar = 1.5*tau21lo - 0.5*tau21hi;
470 
471  Real tau22lo = 0.5 * ( tau22(i , j , k ) + tau22(i , j-1, k ) );
472  Real tau22hi = 0.5 * ( tau22(i , j , k+1) + tau22(i , j-1, k+1) );
473  Real tau22bar = 1.5*tau22lo - 0.5*tau22hi;
474 
475  Real mu_bar = 0.25*( mu_turb(i, j-1, k , EddyDiff::Mom_v) + mu_turb(i, j, k , EddyDiff::Mom_v)
476  + mu_turb(i, j-1, k-1, EddyDiff::Mom_v) + mu_turb(i, j, k-1, EddyDiff::Mom_v) );
477  Real rhoAlpha_bar = 0.25*( rhoAlpha(i, j-1, k ) + rhoAlpha(i, j, k )
478  + rhoAlpha(i, j-1, k-1) + rhoAlpha(i, j, k-1) );
479  Real mu_tot = rhoAlpha_bar + 2.0*mu_bar;
480 
481  tau23(i,j,k) -= met_h_xi*mfx*tau21bar + met_h_eta*mfy*tau22bar;
482  tau23(i,j,k) *= -mu_tot;
483 
484  tau32(i,j,k) *= -mu_tot*met_h_zeta/mfx;
485  });
486  }
487  // Extrapolate tau13 & tau23 to top
488  {
489  Box planexz = tbxxz; planexz.setSmall(2, planexz.bigEnd(2) );
490  tbxxz.growHi(2,-1);
491  ParallelFor(planexz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
492  {
493  Real mfx = mf_ux(i,j,0);
494  Real mfy = mf_uy(i,j,0);
495 
496  Real met_h_xi,met_h_eta,met_h_zeta;
497  met_h_xi = Compute_h_xi_AtEdgeCenterJ (i,j,k,dxInv,z_nd);
498  met_h_eta = Compute_h_eta_AtEdgeCenterJ (i,j,k,dxInv,z_nd);
499  met_h_zeta = Compute_h_zeta_AtEdgeCenterJ(i,j,k,dxInv,z_nd);
500 
501  Real tau11lo = 0.5 * ( tau11(i , j , k-2) + tau11(i-1, j , k-2) );
502  Real tau11hi = 0.5 * ( tau11(i , j , k-1) + tau11(i-1, j , k-1) );
503  Real tau11bar = 1.5*tau11hi - 0.5*tau11lo;
504 
505  Real tau12lo = 0.5 * ( tau12(i , j , k-2) + tau12(i , j+1, k-2) );
506  Real tau12hi = 0.5 * ( tau12(i , j , k-1) + tau12(i , j+1, k-1) );
507  Real tau12bar = 1.5*tau12hi - 0.5*tau12lo;
508 
509  Real mu_bar = 0.25*( mu_turb(i-1, j, k , EddyDiff::Mom_v) + mu_turb(i, j, k , EddyDiff::Mom_v)
510  + mu_turb(i-1, j, k-1, EddyDiff::Mom_v) + mu_turb(i, j, k-1, EddyDiff::Mom_v) );
511  Real rhoAlpha_bar = 0.25*( rhoAlpha(i-1, j, k ) + rhoAlpha(i, j, k )
512  + rhoAlpha(i-1, j, k-1) + rhoAlpha(i, j, k-1) );
513  Real mu_tot = rhoAlpha_bar + 2.0*mu_bar;
514 
515  tau13(i,j,k) -= met_h_xi*mfx*tau11bar + met_h_eta*mfy*tau12bar;
516  tau13(i,j,k) *= -mu_tot;
517 
518  tau31(i,j,k) *= -mu_tot*met_h_zeta/mfy;
519  });
520 
521  Box planeyz = tbxyz; planeyz.setSmall(2, planeyz.bigEnd(2) );
522  tbxyz.growHi(2,-1);
523  ParallelFor(planeyz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
524  {
525  Real mfx = mf_vx(i,j,0);
526  Real mfy = mf_vy(i,j,0);
527 
528  Real met_h_xi,met_h_eta,met_h_zeta;
529  met_h_xi = Compute_h_xi_AtEdgeCenterI (i,j,k,dxInv,z_nd);
530  met_h_eta = Compute_h_eta_AtEdgeCenterI (i,j,k,dxInv,z_nd);
531  met_h_zeta = Compute_h_zeta_AtEdgeCenterI(i,j,k,dxInv,z_nd);
532 
533  Real tau21lo = 0.5 * ( tau21(i , j , k-2) + tau21(i+1, j , k-2) );
534  Real tau21hi = 0.5 * ( tau21(i , j , k-1) + tau21(i+1, j , k-1) );
535  Real tau21bar = 1.5*tau21hi - 0.5*tau21lo;
536 
537  Real tau22lo = 0.5 * ( tau22(i , j , k-2) + tau22(i , j-1, k-2) );
538  Real tau22hi = 0.5 * ( tau22(i , j , k-1) + tau22(i , j-1, k-1) );
539  Real tau22bar = 1.5*tau22hi - 0.5*tau22lo;
540 
541  Real mu_bar = 0.25*( mu_turb(i, j-1, k , EddyDiff::Mom_v) + mu_turb(i, j, k , EddyDiff::Mom_v)
542  + mu_turb(i, j-1, k-1, EddyDiff::Mom_v) + mu_turb(i, j, k-1, EddyDiff::Mom_v) );
543  Real rhoAlpha_bar = 0.25*( rhoAlpha(i, j-1, k ) + rhoAlpha(i, j, k )
544  + rhoAlpha(i, j-1, k-1) + rhoAlpha(i, j, k-1) );
545  Real mu_tot = rhoAlpha_bar + 2.0*mu_bar;
546 
547  tau23(i,j,k) -= met_h_xi*mfx*tau21bar + met_h_eta*mfy*tau22bar;
548  tau23(i,j,k) *= -mu_tot;
549 
550  tau32(i,j,k) *= -mu_tot*met_h_zeta/mfx;
551  });
552  }
553 
554  // Second block: compute 2mu*JT*(S-D)
555  //***********************************************************************************
556  // Fill tau13, tau23 next (valid averaging region)
557  //-----------------------------------------------------------------------------------
558  ParallelFor(tbxxz,tbxyz,
559  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
560  {
561  Real mfx = mf_ux(i,j,0);
562  Real mfy = mf_uy(i,j,0);
563 
564  Real met_h_xi,met_h_eta,met_h_zeta;
565  met_h_xi = Compute_h_xi_AtEdgeCenterJ (i,j,k,dxInv,z_nd);
566  met_h_eta = Compute_h_eta_AtEdgeCenterJ (i,j,k,dxInv,z_nd);
567  met_h_zeta = Compute_h_zeta_AtEdgeCenterJ(i,j,k,dxInv,z_nd);
568 
569  Real tau11bar = 0.25 * ( tau11(i , j , k ) + tau11(i-1, j , k )
570  + tau11(i , j , k-1) + tau11(i-1, j , k-1) );
571  Real tau12bar = 0.25 * ( tau12(i , j , k ) + tau12(i , j+1, k )
572  + tau12(i , j , k-1) + tau12(i , j+1, k-1) );
573 
574  Real mu_bar = 0.25 * ( mu_turb(i-1, j , k , EddyDiff::Mom_v) + mu_turb(i , j , k , EddyDiff::Mom_v)
575  + mu_turb(i-1, j , k-1, EddyDiff::Mom_v) + mu_turb(i , j , k-1, EddyDiff::Mom_v) );
576  Real rhoAlpha_bar = 0.25 * ( rhoAlpha(i-1, j , k ) + rhoAlpha(i , j , k )
577  + rhoAlpha(i-1, j , k-1) + rhoAlpha(i , j , k-1) );
578  Real mu_tot = rhoAlpha_bar + 2.0*mu_bar;
579 
580  tau13(i,j,k) -= met_h_xi*mfx*tau11bar + met_h_eta*mfy*tau12bar;
581  tau13(i,j,k) *= -mu_tot;
582 
583  tau31(i,j,k) *= -mu_tot*met_h_zeta/mfy;
584  },
585  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
586  {
587  Real mfx = mf_vx(i,j,0);
588  Real mfy = mf_vy(i,j,0);
589 
590  Real met_h_xi,met_h_eta,met_h_zeta;
591  met_h_xi = Compute_h_xi_AtEdgeCenterI (i,j,k,dxInv,z_nd);
592  met_h_eta = Compute_h_eta_AtEdgeCenterI (i,j,k,dxInv,z_nd);
593  met_h_zeta = Compute_h_zeta_AtEdgeCenterI(i,j,k,dxInv,z_nd);
594 
595  Real tau21bar = 0.25 * ( tau21(i , j , k ) + tau21(i+1, j , k )
596  + tau21(i , j , k-1) + tau21(i+1, j , k-1) );
597  Real tau22bar = 0.25 * ( tau22(i , j , k ) + tau22(i , j-1, k )
598  + tau22(i , j , k-1) + tau22(i , j-1, k-1) );
599 
600  Real mu_bar = 0.25 * ( 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 = 0.25 * ( 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 + 2.0*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 
609  tau32(i,j,k) *= -mu_tot*met_h_zeta/mfx;
610  });
611 
612  // Fill the remaining components: tau11, tau22, tau12/21
613  //-----------------------------------------------------------------------------------
614  ParallelFor(bxcc,tbxxy,
615  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
616  {
617  Real mfx = mf_mx(i,j,0);
618  Real mfy = mf_my(i,j,0);
619 
620  Real met_h_zeta = detJ(i,j,k);
621 
622  Real mu_tot = rhoAlpha(i,j,k) + 2.0*mu_turb(i, j, k, EddyDiff::Mom_h);
623 
624  tau11(i,j,k) *= -mu_tot*met_h_zeta/mfy;
625  tau22(i,j,k) *= -mu_tot*met_h_zeta/mfx;
626  },
627  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
628  {
629  Real mfx = 0.5 * (mf_ux(i,j,0) + mf_ux(i,j-1,0));
630  Real mfy = 0.5 * (mf_vy(i,j,0) + mf_vy(i-1,j,0));
631 
632  Real met_h_zeta = Compute_h_zeta_AtEdgeCenterK(i,j,k,dxInv,z_nd);
633 
634  Real mu_bar = 0.25*( mu_turb(i-1, j , k, EddyDiff::Mom_h) + mu_turb(i, j , k, EddyDiff::Mom_h)
635  + mu_turb(i-1, j-1, k, EddyDiff::Mom_h) + mu_turb(i, j-1, k, EddyDiff::Mom_h) );
636  Real rhoAlpha_bar = 0.25*( rhoAlpha(i-1, j , k) + rhoAlpha(i, j , k)
637  + rhoAlpha(i-1, j-1, k) + rhoAlpha(i, j-1, k) );
638  Real mu_tot = rhoAlpha_bar + 2.0*mu_bar;
639 
640  tau12(i,j,k) *= -mu_tot*met_h_zeta/mfx;
641  tau21(i,j,k) *= -mu_tot*met_h_zeta/mfy;
642  });
643 }
@ Mom_h
Definition: ERF_IndexDefines.H:170
@ 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: