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 > &tau13, Array4< Real > &tau21, Array4< Real > &tau23, Array4< Real > &tau31, 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)
 
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 > &tau13, Array4< Real > &tau21, Array4< Real > &tau23, Array4< Real > &tau31, 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)
 

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 > &  tau13,
Array4< Real > &  tau21,
Array4< Real > &  tau23,
Array4< Real > &  tau31,
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 
)

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

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 > &  tau13,
Array4< Real > &  tau21,
Array4< Real > &  tau23,
Array4< Real > &  tau31,
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 
)

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
314 {
315  // Handle constant alpha case, in which the provided mu_eff is actually
316  // "alpha" and the viscosity needs to be scaled by rho. This can be further
317  // optimized with if statements below instead of creating a new FAB,
318  // but this is implementation is cleaner.
319  FArrayBox temp;
320  Box gbx = bxcc; // Note: bxcc have been grown in x/y only.
321  gbx.grow(IntVect(0,0,1));
322  temp.resize(gbx,1, The_Async_Arena());
323  Array4<Real> rhoAlpha = temp.array();
324  if (cell_data) {
325  ParallelFor(gbx,
326  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
327  rhoAlpha(i,j,k) = cell_data(i, j, k, Rho_comp) * mu_eff;
328  });
329  } else {
330  ParallelFor(gbx,
331  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
332  rhoAlpha(i,j,k) = mu_eff;
333  });
334  }
335 
336  //***********************************************************************************
337  // NOTE: The first block computes (S-D).
338  // The second block computes 2mu*JT*(S-D)
339  // Boxes are copied here for extrapolations in the second block operations
340  //***********************************************************************************
341  Box bxcc2 = bxcc;
342  bxcc2.grow(IntVect(-1,-1,0));
343 
344  // First block: compute S-D
345  //***********************************************************************************
346  Real OneThird = (1./3.);
347  ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
348  tau11(i,j,k) -= OneThird*er_arr(i,j,k);
349  tau22(i,j,k) -= OneThird*er_arr(i,j,k);
350  tau33(i,j,k) -= OneThird*er_arr(i,j,k);
351  });
352 
353  // Second block: compute 2mu*JT*(S-D)
354  //***********************************************************************************
355  // Fill tau33 first (no linear combination extrapolation)
356  //-----------------------------------------------------------------------------------
357  ParallelFor(bxcc2,
358  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
359  {
360  Real met_h_xi,met_h_eta;
361  met_h_xi = Compute_h_xi_AtCellCenter (i,j,k,dxInv,z_nd);
362  met_h_eta = Compute_h_eta_AtCellCenter (i,j,k,dxInv,z_nd);
363 
364  Real tau31bar = 0.25 * ( tau31(i , j , k ) + tau31(i+1, j , k )
365  + tau31(i , j , k+1) + tau31(i+1, j , k+1) );
366  Real tau32bar = 0.25 * ( tau32(i , j , k ) + tau32(i , j+1, k )
367  + tau32(i , j , k+1) + tau32(i , j+1, k+1) );
368 
369  Real mu_tot = rhoAlpha(i,j,k) + 2.0*mu_turb(i, j, k, EddyDiff::Mom_v);
370 
371  tau33(i,j,k) -= met_h_xi*tau31bar + met_h_eta*tau32bar;
372  tau33(i,j,k) *= -mu_tot;
373  });
374 
375  // Second block: compute 2mu*JT*(S-D)
376  //***********************************************************************************
377  // Fill tau13, tau23 next (linear combination extrapolation)
378  //-----------------------------------------------------------------------------------
379  // Extrapolate tau13 & tau23 to bottom
380  {
381  Box planexz = tbxxz; planexz.setBig(2, planexz.smallEnd(2) );
382  tbxxz.growLo(2,-1);
383  ParallelFor(planexz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
384  {
385  Real met_h_xi,met_h_eta,met_h_zeta;
386  met_h_xi = Compute_h_xi_AtEdgeCenterJ (i,j,k,dxInv,z_nd);
387  met_h_eta = Compute_h_eta_AtEdgeCenterJ (i,j,k,dxInv,z_nd);
388  met_h_zeta = Compute_h_zeta_AtEdgeCenterJ(i,j,k,dxInv,z_nd);
389 
390  Real tau11lo = 0.5 * ( tau11(i , j , k ) + tau11(i-1, j , k ) );
391  Real tau11hi = 0.5 * ( tau11(i , j , k+1) + tau11(i-1, j , k+1) );
392  Real tau11bar = 1.5*tau11lo - 0.5*tau11hi;
393 
394  Real tau12lo = 0.5 * ( tau12(i , j , k ) + tau12(i , j+1, k ) );
395  Real tau12hi = 0.5 * ( tau12(i , j , k+1) + tau12(i , j+1, k+1) );
396  Real tau12bar = 1.5*tau12lo - 0.5*tau12hi;
397 
398  Real mu_bar = 0.25*( mu_turb(i-1, j, k , EddyDiff::Mom_v) + mu_turb(i, j, k , EddyDiff::Mom_v)
399  + mu_turb(i-1, j, k-1, EddyDiff::Mom_v) + mu_turb(i, j, k-1, EddyDiff::Mom_v) );
400  Real rhoAlpha_bar = 0.25*( rhoAlpha(i-1, j, k ) + rhoAlpha(i, j, k )
401  + rhoAlpha(i-1, j, k-1) + rhoAlpha(i, j, k-1) );
402  Real mu_tot = rhoAlpha_bar + 2.0*mu_bar;
403 
404  tau13(i,j,k) -= met_h_xi*tau11bar + met_h_eta*tau12bar;
405  tau13(i,j,k) *= -mu_tot;
406 
407  tau31(i,j,k) *= -mu_tot*met_h_zeta;
408  });
409 
410  Box planeyz = tbxyz; planeyz.setBig(2, planeyz.smallEnd(2) );
411  tbxyz.growLo(2,-1);
412  ParallelFor(planeyz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
413  {
414  Real met_h_xi,met_h_eta,met_h_zeta;
415  met_h_xi = Compute_h_xi_AtEdgeCenterI (i,j,k,dxInv,z_nd);
416  met_h_eta = Compute_h_eta_AtEdgeCenterI (i,j,k,dxInv,z_nd);
417  met_h_zeta = Compute_h_zeta_AtEdgeCenterI(i,j,k,dxInv,z_nd);
418 
419  Real tau21lo = 0.5 * ( tau21(i , j , k ) + tau21(i+1, j , k ) );
420  Real tau21hi = 0.5 * ( tau21(i , j , k+1) + tau21(i+1, j , k+1) );
421  Real tau21bar = 1.5*tau21lo - 0.5*tau21hi;
422 
423  Real tau22lo = 0.5 * ( tau22(i , j , k ) + tau22(i , j-1, k ) );
424  Real tau22hi = 0.5 * ( tau22(i , j , k+1) + tau22(i , j-1, k+1) );
425  Real tau22bar = 1.5*tau22lo - 0.5*tau22hi;
426 
427  Real mu_bar = 0.25*( mu_turb(i, j-1, k , EddyDiff::Mom_v) + mu_turb(i, j, k , EddyDiff::Mom_v)
428  + mu_turb(i, j-1, k-1, EddyDiff::Mom_v) + mu_turb(i, j, k-1, EddyDiff::Mom_v) );
429  Real rhoAlpha_bar = 0.25*( rhoAlpha(i, j-1, k ) + rhoAlpha(i, j, k )
430  + rhoAlpha(i, j-1, k-1) + rhoAlpha(i, j, k-1) );
431  Real mu_tot = rhoAlpha_bar + 2.0*mu_bar;
432 
433  tau23(i,j,k) -= met_h_xi*tau21bar + met_h_eta*tau22bar;
434  tau23(i,j,k) *= -mu_tot;
435 
436  tau32(i,j,k) *= -mu_tot*met_h_zeta;
437  });
438  }
439  // Extrapolate tau13 & tau23 to top
440  {
441  Box planexz = tbxxz; planexz.setSmall(2, planexz.bigEnd(2) );
442  tbxxz.growHi(2,-1);
443  ParallelFor(planexz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
444  {
445  Real met_h_xi,met_h_eta,met_h_zeta;
446  met_h_xi = Compute_h_xi_AtEdgeCenterJ (i,j,k,dxInv,z_nd);
447  met_h_eta = Compute_h_eta_AtEdgeCenterJ (i,j,k,dxInv,z_nd);
448  met_h_zeta = Compute_h_zeta_AtEdgeCenterJ(i,j,k,dxInv,z_nd);
449 
450  Real tau11lo = 0.5 * ( tau11(i , j , k-2) + tau11(i-1, j , k-2) );
451  Real tau11hi = 0.5 * ( tau11(i , j , k-1) + tau11(i-1, j , k-1) );
452  Real tau11bar = 1.5*tau11hi - 0.5*tau11lo;
453 
454  Real tau12lo = 0.5 * ( tau12(i , j , k-2) + tau12(i , j+1, k-2) );
455  Real tau12hi = 0.5 * ( tau12(i , j , k-1) + tau12(i , j+1, k-1) );
456  Real tau12bar = 1.5*tau12hi - 0.5*tau12lo;
457 
458  Real mu_bar = 0.25*( mu_turb(i-1, j, k , EddyDiff::Mom_v) + mu_turb(i, j, k , EddyDiff::Mom_v)
459  + mu_turb(i-1, j, k-1, EddyDiff::Mom_v) + mu_turb(i, j, k-1, EddyDiff::Mom_v) );
460  Real rhoAlpha_bar = 0.25*( rhoAlpha(i-1, j, k ) + rhoAlpha(i, j, k )
461  + rhoAlpha(i-1, j, k-1) + rhoAlpha(i, j, k-1) );
462  Real mu_tot = rhoAlpha_bar + 2.0*mu_bar;
463 
464  tau13(i,j,k) -= met_h_xi*tau11bar + met_h_eta*tau12bar;
465  tau13(i,j,k) *= -mu_tot;
466 
467  tau31(i,j,k) *= -mu_tot*met_h_zeta;
468  });
469 
470  Box planeyz = tbxyz; planeyz.setSmall(2, planeyz.bigEnd(2) );
471  tbxyz.growHi(2,-1);
472  ParallelFor(planeyz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
473  {
474  Real met_h_xi,met_h_eta,met_h_zeta;
475  met_h_xi = Compute_h_xi_AtEdgeCenterI (i,j,k,dxInv,z_nd);
476  met_h_eta = Compute_h_eta_AtEdgeCenterI (i,j,k,dxInv,z_nd);
477  met_h_zeta = Compute_h_zeta_AtEdgeCenterI(i,j,k,dxInv,z_nd);
478 
479  Real tau21lo = 0.5 * ( tau21(i , j , k-2) + tau21(i+1, j , k-2) );
480  Real tau21hi = 0.5 * ( tau21(i , j , k-1) + tau21(i+1, j , k-1) );
481  Real tau21bar = 1.5*tau21hi - 0.5*tau21lo;
482 
483  Real tau22lo = 0.5 * ( tau22(i , j , k-2) + tau22(i , j-1, k-2) );
484  Real tau22hi = 0.5 * ( tau22(i , j , k-1) + tau22(i , j-1, k-1) );
485  Real tau22bar = 1.5*tau22hi - 0.5*tau22lo;
486 
487  Real mu_bar = 0.25*( mu_turb(i, j-1, k , EddyDiff::Mom_v) + mu_turb(i, j, k , EddyDiff::Mom_v)
488  + mu_turb(i, j-1, k-1, EddyDiff::Mom_v) + mu_turb(i, j, k-1, EddyDiff::Mom_v) );
489  Real rhoAlpha_bar = 0.25*( rhoAlpha(i, j-1, k ) + rhoAlpha(i, j, k )
490  + rhoAlpha(i, j-1, k-1) + rhoAlpha(i, j, k-1) );
491  Real mu_tot = rhoAlpha_bar + 2.0*mu_bar;
492 
493  tau23(i,j,k) -= met_h_xi*tau21bar + met_h_eta*tau22bar;
494  tau23(i,j,k) *= -mu_tot;
495 
496  tau32(i,j,k) *= -mu_tot*met_h_zeta;
497  });
498  }
499 
500  // Second block: compute 2mu*JT*(S-D)
501  //***********************************************************************************
502  // Fill tau13, tau23 next (valid averaging region)
503  //-----------------------------------------------------------------------------------
504  ParallelFor(tbxxz,tbxyz,
505  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
506  {
507  Real met_h_xi,met_h_eta,met_h_zeta;
508  met_h_xi = Compute_h_xi_AtEdgeCenterJ (i,j,k,dxInv,z_nd);
509  met_h_eta = Compute_h_eta_AtEdgeCenterJ (i,j,k,dxInv,z_nd);
510  met_h_zeta = Compute_h_zeta_AtEdgeCenterJ(i,j,k,dxInv,z_nd);
511 
512  Real tau11bar = 0.25 * ( tau11(i , j , k ) + tau11(i-1, j , k )
513  + tau11(i , j , k-1) + tau11(i-1, j , k-1) );
514  Real tau12bar = 0.25 * ( tau12(i , j , k ) + tau12(i , j+1, k )
515  + tau12(i , j , k-1) + tau12(i , j+1, k-1) );
516 
517  Real mu_bar = 0.25 * ( mu_turb(i-1, j , k , EddyDiff::Mom_v) + mu_turb(i , j , k , EddyDiff::Mom_v)
518  + mu_turb(i-1, j , k-1, EddyDiff::Mom_v) + mu_turb(i , j , k-1, EddyDiff::Mom_v) );
519  Real rhoAlpha_bar = 0.25 * ( rhoAlpha(i-1, j , k ) + rhoAlpha(i , j , k )
520  + rhoAlpha(i-1, j , k-1) + rhoAlpha(i , j , k-1) );
521  Real mu_tot = rhoAlpha_bar + 2.0*mu_bar;
522 
523  tau13(i,j,k) -= met_h_xi*tau11bar + met_h_eta*tau12bar;
524  tau13(i,j,k) *= -mu_tot;
525 
526  tau31(i,j,k) *= -mu_tot*met_h_zeta;
527  },
528  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
529  {
530  Real met_h_xi,met_h_eta,met_h_zeta;
531  met_h_xi = Compute_h_xi_AtEdgeCenterI (i,j,k,dxInv,z_nd);
532  met_h_eta = Compute_h_eta_AtEdgeCenterI (i,j,k,dxInv,z_nd);
533  met_h_zeta = Compute_h_zeta_AtEdgeCenterI(i,j,k,dxInv,z_nd);
534 
535  Real tau21bar = 0.25 * ( tau21(i , j , k ) + tau21(i+1, j , k )
536  + tau21(i , j , k-1) + tau21(i+1, j , k-1) );
537  Real tau22bar = 0.25 * ( tau22(i , j , k ) + tau22(i , j-1, k )
538  + tau22(i , j , k-1) + tau22(i , j-1, k-1) );
539 
540  Real mu_bar = 0.25 * ( mu_turb(i , j-1, k , EddyDiff::Mom_v) + mu_turb(i , j , k , EddyDiff::Mom_v)
541  + mu_turb(i , j-1, k-1, EddyDiff::Mom_v) + mu_turb(i , j , k-1, EddyDiff::Mom_v) );
542  Real rhoAlpha_bar = 0.25 * ( rhoAlpha(i , j-1, k ) + rhoAlpha(i , j , k )
543  + rhoAlpha(i , j-1, k-1) + rhoAlpha(i , j , k-1) );
544  Real mu_tot = rhoAlpha_bar + 2.0*mu_bar;
545 
546  tau23(i,j,k) -= met_h_xi*tau21bar + met_h_eta*tau22bar;
547  tau23(i,j,k) *= -mu_tot;
548 
549  tau32(i,j,k) *= -mu_tot*met_h_zeta;
550  });
551 
552  // Fill the remaining components: tau11, tau22, tau12/21
553  //-----------------------------------------------------------------------------------
554  ParallelFor(bxcc,tbxxy,
555  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
556  {
557  Real met_h_zeta = detJ(i,j,k);
558 
559  Real mu_tot = rhoAlpha(i,j,k) + 2.0*mu_turb(i, j, k, EddyDiff::Mom_h);
560 
561  tau11(i,j,k) *= -mu_tot*met_h_zeta;
562  tau22(i,j,k) *= -mu_tot*met_h_zeta;
563  },
564  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
565  {
566  Real met_h_zeta = Compute_h_zeta_AtEdgeCenterK(i,j,k,dxInv,z_nd);
567 
568  Real mu_bar = 0.25*( mu_turb(i-1, j , k, EddyDiff::Mom_h) + mu_turb(i, j , k, EddyDiff::Mom_h)
569  + mu_turb(i-1, j-1, k, EddyDiff::Mom_h) + mu_turb(i, j-1, k, EddyDiff::Mom_h) );
570  Real rhoAlpha_bar = 0.25*( rhoAlpha(i-1, j , k) + rhoAlpha(i, j , k)
571  + rhoAlpha(i-1, j-1, k) + rhoAlpha(i, j-1, k) );
572  Real mu_tot = rhoAlpha_bar + 2.0*mu_bar;
573 
574  tau12(i,j,k) *= -mu_tot*met_h_zeta;
575  tau21(i,j,k) *= -mu_tot*met_h_zeta;
576  });
577 }
EddyDiff::Mom_h
@ Mom_h
Definition: ERF_IndexDefines.H:151
EddyDiff::Mom_v
@ Mom_v
Definition: ERF_IndexDefines.H:156

Referenced by erf_make_tau_terms().

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