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 <ERF_TerrainImplicitMetric.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, const bool implicit_metric)
 
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, const bool implicit_metric)
 

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,
const bool  implicit_metric 
)

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

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,
const bool  implicit_metric 
)

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
[in]implicit_metricerf.implicit_terrain_metric: tau13i/tau23i also hold the compact terrain-metric term on the interior faces (see ERF_TerrainImplicitMetric.H)

NOTE: The zeta-face stresses tau13/tau23 are the terrain-normal combination tau_i3 - h_xi*tau_i1 - h_eta*tau_i2. tau_i3 is a vertical stress and takes K_v (EddyDiff::Mom_v); tau_i1 and tau_i2 are horizontal stresses and take K_h (EddyDiff::Mom_h), as the h*Fx terms of the scalar fluxes do. Applying K_v to the projected terms makes the operator non-symmetric and, when K_h*h^2 > 2*K_v, anti-diffusive. The projected stresses K_h*S are formed at the cells and xy edges and averaged to the zeta edges. On a uniform slope with uniform dz this is the transpose of the metric term in S11/S22, so the strain part of the operator dissipates energy for any K_h and K_v (an edge-averaged K_h times averaged strains does not, where K_h varies from cell to cell). tau13i/tau23i (the part the implicit solve takes) are pure K_v terms.

428 {
429  // NOTE: mu_eff includes factor of 2
430 
431  // Handle constant alpha case, in which the provided mu_eff is actually
432  // "alpha" and the viscosity needs to be scaled by rho. This can be further
433  // optimized with if statements below instead of creating a new FAB,
434  // but this is implementation is cleaner.
435  FArrayBox temp;
436  Box gbx = bxcc; // Note: bxcc have been grown in x/y only.
437  gbx.grow(IntVect(0,0,1));
438  temp.resize(gbx,1, The_Async_Arena());
439  Array4<Real> rhoAlpha = temp.array();
440 
441  if (cell_data)
442  // constant alpha (stored in mu_eff)
443  {
444  ParallelFor(gbx,
445  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
446  rhoAlpha(i,j,k) = cell_data(i, j, k, Rho_comp) * mu_eff;
447  });
448  }
449  else
450  // constant mu_eff
451  {
452  ParallelFor(gbx,
453  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
454  rhoAlpha(i,j,k) = mu_eff;
455  });
456  }
457 
458  //***********************************************************************************
459  // NOTE: The first block computes (S-D).
460  // The second block computes JT*2K*(S-D)
461  // Boxes are copied here for extrapolations in the second block operations
462  //***********************************************************************************
463  Box bxcc2 = bxcc;
464  bxcc2.grow(IntVect(-1,-1,0));
465 
466  // First block: compute S-D
467  //***********************************************************************************
468  Real OneThird = (one/three);
469  ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
470  if (tau33i) tau33i(i,j,k) = tau33(i,j,k);
471 
472  tau11(i,j,k) -= OneThird*er_arr(i,j,k);
473  tau22(i,j,k) -= OneThird*er_arr(i,j,k);
474  tau33(i,j,k) -= OneThird*er_arr(i,j,k);
475  });
476 
477  // Second block: compute JT*2K*(S-D)
478  //***********************************************************************************
479  // Fill tau33 first (no linear combination extrapolation)
480  //-----------------------------------------------------------------------------------
481  ParallelFor(bxcc2,
482  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
483  {
484  Real mfx = mf_mx(i,j,0);
485  Real mfy = mf_my(i,j,0);
486 
487  Real met_h_xi,met_h_eta;
488  met_h_xi = Compute_h_xi_AtCellCenter (i,j,k,dxInv,z_nd);
489  met_h_eta = Compute_h_eta_AtCellCenter (i,j,k,dxInv,z_nd);
490 
491  Real tau31bar = fourth * ( tau31(i , j , k ) + tau31(i+1, j , k )
492  + tau31(i , j , k+1) + tau31(i+1, j , k+1) );
493  Real tau32bar = fourth * ( tau32(i , j , k ) + tau32(i , j+1, k )
494  + tau32(i , j , k+1) + tau32(i , j+1, k+1) );
495 
496  Real mu_tot = rhoAlpha(i,j,k) + two*mu_turb(i, j, k, EddyDiff::Mom_v);
497 
498  tau33(i,j,k) -= met_h_xi*mfx*tau31bar + met_h_eta*mfy*tau32bar;
499  tau33(i,j,k) *= -mu_tot;
500 
501  if (tau33i) { tau33i(i,j,k) *= -mu_tot; }
502  });
503 
504  // The zeta-face stresses project the horizontal stresses K_h*S11, K_h*S12, K_h*S21 and
505  // K_h*S22 (see the NOTE above). Form them where they live -- at the cells and the xy
506  // edges, with the coefficients the tau11/tau22/tau12 kernels below apply -- and average
507  // the stresses to the zeta edges, as the scalar fluxes average K_h*grad. On a uniform
508  // slope with uniform dz that is the transpose of the metric term in S11/S22, which keeps
509  // the strain part of the operator dissipative for any K_h and K_v (see the NOTE above).
510  //-----------------------------------------------------------------------------------
511  // kh12 is read at (i,j+1) over tbxxz and kh21 at (i+1,j) over tbxyz. On a tile that is not
512  // the last in x or y the nodal tile box tbxxy stops one node short of those reads, so form
513  // them on the nodes the reads need; the strains are valid there (they are computed on the
514  // halo-grown boxes).
515  Box khxy = tbxxy;
516  khxy.setSmall(0, amrex::min(tbxxz.smallEnd(0), tbxyz.smallEnd(0)));
517  khxy.setSmall(1, amrex::min(tbxxz.smallEnd(1), tbxyz.smallEnd(1)));
518  khxy.setBig (0, amrex::max(tbxxz.bigEnd(0) , tbxyz.bigEnd(0)+1));
519  khxy.setBig (1, amrex::max(tbxxz.bigEnd(1)+1, tbxyz.bigEnd(1) ));
520  FArrayBox kh11_fab(bxcc,1,The_Async_Arena()), kh22_fab(bxcc,1,The_Async_Arena());
521  FArrayBox kh12_fab(khxy,1,The_Async_Arena()), kh21_fab(khxy,1,The_Async_Arena());
522  const Array4<Real> kh11 = kh11_fab.array();
523  const Array4<Real> kh22 = kh22_fab.array();
524  const Array4<Real> kh12 = kh12_fab.array();
525  const Array4<Real> kh21 = kh21_fab.array();
526  ParallelFor(bxcc, khxy,
527  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
528  {
529  Real mu_h_tot = rhoAlpha(i,j,k) + two*mu_turb(i, j, k, EddyDiff::Mom_h);
530  kh11(i,j,k) = mu_h_tot*tau11(i,j,k);
531  kh22(i,j,k) = mu_h_tot*tau22(i,j,k);
532  },
533  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
534  {
535  Real mu_bar = fourth * ( mu_turb(i-1, j , k, EddyDiff::Mom_h) + mu_turb(i, j , k, EddyDiff::Mom_h)
536  + mu_turb(i-1, j-1, k, EddyDiff::Mom_h) + mu_turb(i, j-1, k, EddyDiff::Mom_h) );
537  Real rhoAlpha_bar = fourth * ( rhoAlpha(i-1, j , k) + rhoAlpha(i, j , k)
538  + rhoAlpha(i-1, j-1, k) + rhoAlpha(i, j-1, k) );
539  Real mu_h_tot = rhoAlpha_bar + two*mu_bar;
540  kh12(i,j,k) = mu_h_tot*tau12(i,j,k);
541  kh21(i,j,k) = mu_h_tot*tau21(i,j,k);
542  });
543 
544  // Second block: compute JT*2K*(S-D)
545  //***********************************************************************************
546  // Fill tau13, tau23 next (linear combination extrapolation)
547  //-----------------------------------------------------------------------------------
548  // Extrapolate tau13 & tau23 to bottom
549  {
550  Box planexz = tbxxz; planexz.setBig(2, planexz.smallEnd(2) );
551  tbxxz.growLo(2,-1);
552  Box planeyz = tbxyz; planeyz.setBig(2, planeyz.smallEnd(2) );
553  tbxyz.growLo(2,-1);
554  ParallelFor(planexz,planeyz,
555  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
556  {
557  Real mfx = mf_ux(i,j,0);
558  Real mfy = mf_uy(i,j,0);
559 
560  Real met_h_xi,met_h_eta,met_h_zeta;
561  met_h_xi = Compute_h_xi_AtEdgeCenterJ (i,j,k,dxInv,z_nd);
562  met_h_eta = Compute_h_eta_AtEdgeCenterJ (i,j,k,dxInv,z_nd);
563  met_h_zeta = Compute_h_zeta_AtEdgeCenterJ(i,j,k,dxInv,z_nd);
564 
565  Real tau11lo = myhalf * ( kh11(i , j , k ) + kh11(i-1, j , k ) );
566  Real tau11hi = myhalf * ( kh11(i , j , k+1) + kh11(i-1, j , k+1) );
567  Real tau11bar = Real(1.5)*tau11lo - myhalf*tau11hi;
568 
569  Real tau12lo = myhalf * ( kh12(i , j , k ) + kh12(i , j+1, k ) );
570  Real tau12hi = myhalf * ( kh12(i , j , k+1) + kh12(i , j+1, k+1) );
571  Real tau12bar = Real(1.5)*tau12lo - myhalf*tau12hi;
572 
573  Real mu_bar = fourth * ( mu_turb(i-1, j, k , EddyDiff::Mom_v) + mu_turb(i, j, k , EddyDiff::Mom_v)
574  + mu_turb(i-1, j, k-1, EddyDiff::Mom_v) + mu_turb(i, j, k-1, EddyDiff::Mom_v) );
575  Real rhoAlpha_bar = fourth * ( rhoAlpha(i-1, j, k ) + rhoAlpha(i, j, k )
576  + rhoAlpha(i-1, j, k-1) + rhoAlpha(i, j, k-1) );
577  Real mu_tot = rhoAlpha_bar + two*mu_bar;
578 
579  // K_v on S13; the projected horizontal stresses carry K_h (see the NOTE above)
580  tau13(i,j,k) *= -mu_tot;
581  tau13(i,j,k) += met_h_xi*mfx*tau11bar + met_h_eta*mfy*tau12bar;
582  if (tau13i) { tau13i(i,j,k) *= -mu_tot; }
583 
584  tau31(i,j,k) *= -mu_tot*met_h_zeta/mfy;
585  },
586  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
587  {
588  Real mfx = mf_vx(i,j,0);
589  Real mfy = mf_vy(i,j,0);
590 
591  Real met_h_xi,met_h_eta,met_h_zeta;
592  met_h_xi = Compute_h_xi_AtEdgeCenterI (i,j,k,dxInv,z_nd);
593  met_h_eta = Compute_h_eta_AtEdgeCenterI (i,j,k,dxInv,z_nd);
594  met_h_zeta = Compute_h_zeta_AtEdgeCenterI(i,j,k,dxInv,z_nd);
595 
596  Real tau21lo = myhalf * ( kh21(i , j , k ) + kh21(i+1, j , k ) );
597  Real tau21hi = myhalf * ( kh21(i , j , k+1) + kh21(i+1, j , k+1) );
598  Real tau21bar = Real(1.5)*tau21lo - myhalf*tau21hi;
599 
600  Real tau22lo = myhalf * ( kh22(i , j , k ) + kh22(i , j-1, k ) );
601  Real tau22hi = myhalf * ( kh22(i , j , k+1) + kh22(i , j-1, k+1) );
602  Real tau22bar = Real(1.5)*tau22lo - myhalf*tau22hi;
603 
604  Real mu_bar = fourth * ( mu_turb(i, j-1, k , EddyDiff::Mom_v) + mu_turb(i, j, k , EddyDiff::Mom_v)
605  + mu_turb(i, j-1, k-1, EddyDiff::Mom_v) + mu_turb(i, j, k-1, EddyDiff::Mom_v) );
606  Real rhoAlpha_bar = fourth * ( rhoAlpha(i, j-1, k ) + rhoAlpha(i, j, k )
607  + rhoAlpha(i, j-1, k-1) + rhoAlpha(i, j, k-1) );
608  Real mu_tot = rhoAlpha_bar + two*mu_bar;
609 
610  // K_v on S23; the projected horizontal stresses carry K_h (see the NOTE above)
611  tau23(i,j,k) *= -mu_tot;
612  tau23(i,j,k) += met_h_xi*mfx*tau21bar + met_h_eta*mfy*tau22bar;
613  if (tau23i) { tau23i(i,j,k) *= -mu_tot; }
614 
615  tau32(i,j,k) *= -mu_tot*met_h_zeta/mfx;
616  });
617  }
618  // Extrapolate tau13 & tau23 to top
619  {
620  Box planexz = tbxxz; planexz.setSmall(2, planexz.bigEnd(2) );
621  tbxxz.growHi(2,-1);
622  Box planeyz = tbxyz; planeyz.setSmall(2, planeyz.bigEnd(2) );
623  tbxyz.growHi(2,-1);
624  ParallelFor(planexz,planeyz,
625  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
626  {
627  Real mfx = mf_ux(i,j,0);
628  Real mfy = mf_uy(i,j,0);
629 
630  Real met_h_xi,met_h_eta,met_h_zeta;
631  met_h_xi = Compute_h_xi_AtEdgeCenterJ (i,j,k,dxInv,z_nd);
632  met_h_eta = Compute_h_eta_AtEdgeCenterJ (i,j,k,dxInv,z_nd);
633  met_h_zeta = Compute_h_zeta_AtEdgeCenterJ(i,j,k,dxInv,z_nd);
634 
635  Real tau11lo = myhalf * ( kh11(i , j , k-2) + kh11(i-1, j , k-2) );
636  Real tau11hi = myhalf * ( kh11(i , j , k-1) + kh11(i-1, j , k-1) );
637  Real tau11bar = Real(1.5)*tau11hi - myhalf*tau11lo;
638 
639  Real tau12lo = myhalf * ( kh12(i , j , k-2) + kh12(i , j+1, k-2) );
640  Real tau12hi = myhalf * ( kh12(i , j , k-1) + kh12(i , j+1, k-1) );
641  Real tau12bar = Real(1.5)*tau12hi - myhalf*tau12lo;
642 
643  Real mu_bar = fourth * ( mu_turb(i-1, j, k , EddyDiff::Mom_v) + mu_turb(i, j, k , EddyDiff::Mom_v)
644  + mu_turb(i-1, j, k-1, EddyDiff::Mom_v) + mu_turb(i, j, k-1, EddyDiff::Mom_v) );
645  Real rhoAlpha_bar = fourth * ( rhoAlpha(i-1, j, k ) + rhoAlpha(i, j, k )
646  + rhoAlpha(i-1, j, k-1) + rhoAlpha(i, j, k-1) );
647  Real mu_tot = rhoAlpha_bar + two*mu_bar;
648 
649  // K_v on S13; the projected horizontal stresses carry K_h (see the NOTE above)
650  tau13(i,j,k) *= -mu_tot;
651  tau13(i,j,k) += met_h_xi*mfx*tau11bar + met_h_eta*mfy*tau12bar;
652  if (tau13i) { tau13i(i,j,k) *= -mu_tot; }
653 
654  tau31(i,j,k) *= -mu_tot*met_h_zeta/mfy;
655  },
656  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
657  {
658  Real mfx = mf_vx(i,j,0);
659  Real mfy = mf_vy(i,j,0);
660 
661  Real met_h_xi,met_h_eta,met_h_zeta;
662  met_h_xi = Compute_h_xi_AtEdgeCenterI (i,j,k,dxInv,z_nd);
663  met_h_eta = Compute_h_eta_AtEdgeCenterI (i,j,k,dxInv,z_nd);
664  met_h_zeta = Compute_h_zeta_AtEdgeCenterI(i,j,k,dxInv,z_nd);
665 
666  Real tau21lo = myhalf * ( kh21(i , j , k-2) + kh21(i+1, j , k-2) );
667  Real tau21hi = myhalf * ( kh21(i , j , k-1) + kh21(i+1, j , k-1) );
668  Real tau21bar = Real(1.5)*tau21hi - myhalf*tau21lo;
669 
670  Real tau22lo = myhalf * ( kh22(i , j , k-2) + kh22(i , j-1, k-2) );
671  Real tau22hi = myhalf * ( kh22(i , j , k-1) + kh22(i , j-1, k-1) );
672  Real tau22bar = Real(1.5)*tau22hi - myhalf*tau22lo;
673 
674  Real mu_bar = fourth * ( mu_turb(i, j-1, k , EddyDiff::Mom_v) + mu_turb(i, j, k , EddyDiff::Mom_v)
675  + mu_turb(i, j-1, k-1, EddyDiff::Mom_v) + mu_turb(i, j, k-1, EddyDiff::Mom_v) );
676  Real rhoAlpha_bar = fourth * ( rhoAlpha(i, j-1, k ) + rhoAlpha(i, j, k )
677  + rhoAlpha(i, j-1, k-1) + rhoAlpha(i, j, k-1) );
678  Real mu_tot = rhoAlpha_bar + two*mu_bar;
679 
680  // K_v on S23; the projected horizontal stresses carry K_h (see the NOTE above)
681  tau23(i,j,k) *= -mu_tot;
682  tau23(i,j,k) += met_h_xi*mfx*tau21bar + met_h_eta*mfy*tau22bar;
683  if (tau23i) { tau23i(i,j,k) *= -mu_tot; }
684 
685  tau32(i,j,k) *= -mu_tot*met_h_zeta/mfx;
686  });
687  }
688 
689  // Second block: compute JT*2K*(S-D)
690  //***********************************************************************************
691  // Fill tau13, tau23 next (valid averaging region)
692  //-----------------------------------------------------------------------------------
693  ParallelFor(tbxxz,tbxyz,
694  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
695  {
696  Real mfx = mf_ux(i,j,0);
697  Real mfy = mf_uy(i,j,0);
698 
699  Real met_h_xi,met_h_eta,met_h_zeta;
700  met_h_xi = Compute_h_xi_AtEdgeCenterJ (i,j,k,dxInv,z_nd);
701  met_h_eta = Compute_h_eta_AtEdgeCenterJ (i,j,k,dxInv,z_nd);
702  met_h_zeta = Compute_h_zeta_AtEdgeCenterJ(i,j,k,dxInv,z_nd);
703 
704  Real tau11bar = fourth * ( kh11(i , j , k ) + kh11(i-1, j , k )
705  + kh11(i , j , k-1) + kh11(i-1, j , k-1) );
706  Real tau12bar = fourth * ( kh12(i , j , k ) + kh12(i , j+1, k )
707  + kh12(i , j , k-1) + kh12(i , j+1, k-1) );
708 
709  Real mu_bar = fourth * ( mu_turb(i-1, j , k , EddyDiff::Mom_v) + mu_turb(i , j , k , EddyDiff::Mom_v)
710  + mu_turb(i-1, j , k-1, EddyDiff::Mom_v) + mu_turb(i , j , k-1, EddyDiff::Mom_v) );
711  Real rhoAlpha_bar = fourth * ( rhoAlpha(i-1, j , k ) + rhoAlpha(i , j , k )
712  + rhoAlpha(i-1, j , k-1) + rhoAlpha(i , j , k-1) );
713  Real mu_tot = rhoAlpha_bar + two*mu_bar;
714 
715  // K_v on S13; the projected horizontal stresses carry K_h (see the NOTE above)
716  tau13(i,j,k) *= -mu_tot;
717  tau13(i,j,k) += met_h_xi*mfx*tau11bar + met_h_eta*mfy*tau12bar;
718  if (tau13i) {
719  // erf.implicit_terrain_metric: the implicit solve also takes the compact
720  // K_h*M*d/dz metric term on these interior faces, with K_h averaged to the edge
721  // as getRhoAlphaForFaces averages it there, so tau13i holds it too
722  Real metric_coef = zero;
723  if (implicit_metric) {
724  Real mu_h_bar = fourth * ( mu_turb(i-1, j , k , EddyDiff::Mom_h) + mu_turb(i , j , k , EddyDiff::Mom_h)
725  + mu_turb(i-1, j , k-1, EddyDiff::Mom_h) + mu_turb(i , j , k-1, EddyDiff::Mom_h) );
726  metric_coef = (rhoAlpha_bar + two*mu_h_bar)
727  * TerrainMetricMomFactor(0, met_h_xi, met_h_eta, mfx, mfy);
728  }
729  tau13i(i,j,k) *= -(mu_tot + metric_coef);
730  }
731 
732  tau31(i,j,k) *= -mu_tot*met_h_zeta/mfy;
733  },
734  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
735  {
736  Real mfx = mf_vx(i,j,0);
737  Real mfy = mf_vy(i,j,0);
738 
739  Real met_h_xi,met_h_eta,met_h_zeta;
740  met_h_xi = Compute_h_xi_AtEdgeCenterI (i,j,k,dxInv,z_nd);
741  met_h_eta = Compute_h_eta_AtEdgeCenterI (i,j,k,dxInv,z_nd);
742  met_h_zeta = Compute_h_zeta_AtEdgeCenterI(i,j,k,dxInv,z_nd);
743 
744  Real tau21bar = fourth * ( kh21(i , j , k ) + kh21(i+1, j , k )
745  + kh21(i , j , k-1) + kh21(i+1, j , k-1) );
746  Real tau22bar = fourth * ( kh22(i , j , k ) + kh22(i , j-1, k )
747  + kh22(i , j , k-1) + kh22(i , j-1, k-1) );
748 
749  Real mu_bar = fourth * ( mu_turb(i , j-1, k , EddyDiff::Mom_v) + mu_turb(i , j , k , EddyDiff::Mom_v)
750  + mu_turb(i , j-1, k-1, EddyDiff::Mom_v) + mu_turb(i , j , k-1, EddyDiff::Mom_v) );
751  Real rhoAlpha_bar = fourth * ( rhoAlpha(i , j-1, k ) + rhoAlpha(i , j , k )
752  + rhoAlpha(i , j-1, k-1) + rhoAlpha(i , j , k-1) );
753  Real mu_tot = rhoAlpha_bar + two*mu_bar;
754 
755  // K_v on S23; the projected horizontal stresses carry K_h (see the NOTE above)
756  tau23(i,j,k) *= -mu_tot;
757  tau23(i,j,k) += met_h_xi*mfx*tau21bar + met_h_eta*mfy*tau22bar;
758  if (tau23i) {
759  // erf.implicit_terrain_metric: the implicit solve also takes the compact
760  // K_h*M*d/dz metric term on these interior faces, with K_h averaged to the edge
761  // as getRhoAlphaForFaces averages it there, so tau23i holds it too
762  Real metric_coef = zero;
763  if (implicit_metric) {
764  Real mu_h_bar = fourth * ( mu_turb(i , j-1, k , EddyDiff::Mom_h) + mu_turb(i , j , k , EddyDiff::Mom_h)
765  + mu_turb(i , j-1, k-1, EddyDiff::Mom_h) + mu_turb(i , j , k-1, EddyDiff::Mom_h) );
766  metric_coef = (rhoAlpha_bar + two*mu_h_bar)
767  * TerrainMetricMomFactor(1, met_h_xi, met_h_eta, mfx, mfy);
768  }
769  tau23i(i,j,k) *= -(mu_tot + metric_coef);
770  }
771 
772  tau32(i,j,k) *= -mu_tot*met_h_zeta/mfx;
773  });
774 
775  // Fill the remaining components: tau11, tau22, tau12/21
776  //-----------------------------------------------------------------------------------
777  ParallelFor(bxcc,tbxxy,
778  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
779  {
780  Real mfx = mf_mx(i,j,0);
781  Real mfy = mf_my(i,j,0);
782 
783  Real met_h_zeta = detJ(i,j,k);
784 
785  Real mu_tot = rhoAlpha(i,j,k) + two*mu_turb(i, j, k, EddyDiff::Mom_h);
786 
787  tau11(i,j,k) *= -mu_tot*met_h_zeta/mfy;
788  tau22(i,j,k) *= -mu_tot*met_h_zeta/mfx;
789  },
790  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
791  {
792  Real mfx = myhalf * (mf_ux(i,j,0) + mf_ux(i,j-1,0));
793  Real mfy = myhalf * (mf_vy(i,j,0) + mf_vy(i-1,j,0));
794 
795  Real met_h_zeta = Compute_h_zeta_AtEdgeCenterK(i,j,k,dxInv,z_nd);
796 
797  Real mu_bar = fourth * ( mu_turb(i-1, j , k, EddyDiff::Mom_h) + mu_turb(i, j , k, EddyDiff::Mom_h)
798  + mu_turb(i-1, j-1, k, EddyDiff::Mom_h) + mu_turb(i, j-1, k, EddyDiff::Mom_h) );
799  Real rhoAlpha_bar = fourth * ( rhoAlpha(i-1, j , k) + rhoAlpha(i, j , k)
800  + rhoAlpha(i-1, j-1, k) + rhoAlpha(i, j-1, k) );
801  Real mu_tot = rhoAlpha_bar + two*mu_bar;
802 
803  tau12(i,j,k) *= -mu_tot*met_h_zeta/mfx;
804  tau21(i,j,k) *= -mu_tot*met_h_zeta/mfy;
805  });
806 }
two
constexpr amrex::Real two
Definition: ERF_NumericalConstants.H:38
EddyDiff::Mom_h
@ Mom_h
Definition: ERF_IndexDefines.H:247
EddyDiff::Mom_v
@ Mom_v
Definition: ERF_IndexDefines.H:252

Referenced by erf_make_tau_terms().

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