ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
ERF_MakeTauTerms.cpp File Reference
#include "AMReX_ArrayLim.H"
#include "AMReX_BCRec.H"
#include "AMReX_GpuContainers.H"
#include "ERF_TI_slow_headers.H"
#include "ERF_EOS.H"
#include "ERF_Utils.H"
Include dependency graph for ERF_MakeTauTerms.cpp:

Functions

void erf_make_tau_terms (int level, int nrk, const Vector< BCRec > &domain_bcs_type_h, const MultiFab &z_phys_nd, Vector< MultiFab > &S_data, const MultiFab &xvel, const MultiFab &yvel, const MultiFab &zvel, Vector< std::unique_ptr< MultiFab >> &Tau_lev, MultiFab *SmnSmn, MultiFab *eddyDiffs, const Geometry geom, const SolverChoice &solverChoice, std::unique_ptr< SurfaceLayer > &, Gpu::DeviceVector< Real > &stretched_dz_d, const MultiFab &detJ, Vector< std::unique_ptr< MultiFab >> &mapfac)
 

Function Documentation

◆ erf_make_tau_terms()

void erf_make_tau_terms ( int  level,
int  nrk,
const Vector< BCRec > &  domain_bcs_type_h,
const MultiFab &  z_phys_nd,
Vector< MultiFab > &  S_data,
const MultiFab &  xvel,
const MultiFab &  yvel,
const MultiFab &  zvel,
Vector< std::unique_ptr< MultiFab >> &  Tau_lev,
MultiFab *  SmnSmn,
MultiFab *  eddyDiffs,
const Geometry  geom,
const SolverChoice solverChoice,
std::unique_ptr< SurfaceLayer > &  ,
Gpu::DeviceVector< Real > &  stretched_dz_d,
const MultiFab &  detJ,
Vector< std::unique_ptr< MultiFab >> &  mapfac 
)
28 {
29  BL_PROFILE_REGION("erf_make_tau_terms()");
30 
31  const BCRec* bc_ptr_h = domain_bcs_type_h.data();
32 
33  DiffChoice dc = solverChoice.diffChoice;
34  TurbChoice tc = solverChoice.turbChoice[level];
35 
36  const bool l_use_terrain_fitted_coords = (solverChoice.mesh_type != MeshType::ConstantDz);
37  const bool l_moving_terrain = (solverChoice.terrain_type == TerrainType::MovingFittedMesh);
38  if (l_moving_terrain) AMREX_ALWAYS_ASSERT (l_use_terrain_fitted_coords);
39 
40 
41  const bool l_use_diff = ( (dc.molec_diff_type != MolecDiffType::None) || tc.use_kturb );
42  const bool l_use_constAlpha = ( dc.molec_diff_type == MolecDiffType::ConstantAlpha );
43  const bool l_use_turb = ( tc.les_type == LESType::Smagorinsky ||
44  tc.les_type == LESType::Deardorff ||
45  tc.rans_type == RANSType::kEqn ||
46  tc.pbl_type == PBLType::MYJ ||
47  tc.pbl_type == PBLType::MYNN25 ||
48  tc.pbl_type == PBLType::MYNNEDMF ||
49  tc.pbl_type == PBLType::YSU ||
50  tc.pbl_type == PBLType::MRF);
51 
52  const bool need_SmnSmn = (tc.les_type == LESType::Deardorff ||
53  tc.rans_type == RANSType::kEqn);
54 
55  const Box& domain = geom.Domain();
56  const int domlo_z = domain.smallEnd(2);
57  const int domhi_z = domain.bigEnd(2);
58 
59  const GpuArray<Real, AMREX_SPACEDIM> dxInv = geom.InvCellSizeArray();
60 
61  // *****************************************************************************
62  // Pre-computed quantities
63  // *****************************************************************************
64  const BoxArray& ba = S_data[IntVars::cons].boxArray();
65  const DistributionMapping& dm = S_data[IntVars::cons].DistributionMap();
66 
67  std::unique_ptr<MultiFab> expr;
68 
69  if (l_use_diff) {
70  expr = std::make_unique<MultiFab>(ba, dm, 1, IntVect(1,1,1));
71 
72  // if using constant alpha (mu = rho * alpha), then first divide by the
73  // reference density -- mu_eff will be scaled by the instantaneous
74  // local density later when ComputeStress*Visc_*() is called
75  Real mu_eff = (l_use_constAlpha) ? 2.0 * dc.dynamic_viscosity / dc.rho0_trans
76  : 2.0 * dc.dynamic_viscosity;
77 
78  auto dz_ptr = stretched_dz_d.data();
79 
80 #ifdef _OPENMP
81 #pragma omp parallel if (Gpu::notInLaunchRegion())
82 #endif
83  for ( MFIter mfi(S_data[IntVars::cons],TileNoZ()); mfi.isValid(); ++mfi)
84  {
85  const Box& bx = mfi.tilebox();
86  const Box& valid_bx = mfi.validbox();
87 
88  // Velocities
89  const Array4<const Real> & u = xvel.array(mfi);
90  const Array4<const Real> & v = yvel.array(mfi);
91  const Array4<const Real> & w = zvel.array(mfi);
92 
93  // Map factors
94  const Array4<const Real>& mf_mx = mapfac[MapFacType::m_x]->const_array(mfi);
95  const Array4<const Real>& mf_ux = mapfac[MapFacType::u_x]->const_array(mfi);
96  const Array4<const Real>& mf_vx = mapfac[MapFacType::v_x]->const_array(mfi);
97  const Array4<const Real>& mf_my = mapfac[MapFacType::m_y]->const_array(mfi);
98  const Array4<const Real>& mf_uy = mapfac[MapFacType::u_y]->const_array(mfi);
99  const Array4<const Real>& mf_vy = mapfac[MapFacType::v_y]->const_array(mfi);
100 
101  // Eddy viscosity
102  const Array4<Real const>& mu_turb = l_use_turb ? eddyDiffs->const_array(mfi) :
103  Array4<const Real>{};
104  const Array4<Real const>& cell_data = l_use_constAlpha ? S_data[IntVars::cons].const_array(mfi) :
105  Array4<const Real>{};
106 
107  // Terrain metrics
108  const Array4<const Real>& z_nd = z_phys_nd.const_array(mfi);
109  const Array4<const Real>& detJ_arr = detJ.const_array(mfi);
110 
111  //-------------------------------------------------------------------------------
112  // NOTE: Tile boxes with terrain are not intuitive. The linear combination of
113  // stress terms requires care. Create a tile box that intersects the
114  // valid box, then grow the box in x/y. Compute the strain on the local
115  // FAB over this grown tile box. Compute the stress over the tile box,
116  // except tau_ii which still needs the halo cells. Finally, write from
117  // the local FAB to the Tau MF but only on the tile box.
118  //-------------------------------------------------------------------------------
119 
120  //-------------------------------------------------------------------------------
121  // TODO: Avoid recomputing strain on the first RK stage. One could populate
122  // the FABs with tau_ij, compute stress, and then write to tau_ij. The
123  // problem with this approach is you will over-write the needed halo layer
124  // needed by subsequent tile boxes (particularly S_ii becomes Tau_ii).
125  //-------------------------------------------------------------------------------
126 
127  // Strain/Stress tile boxes
128  Box bxcc = mfi.tilebox();
129  Box tbxxy = mfi.tilebox(IntVect(1,1,0));
130  Box tbxxz = mfi.tilebox(IntVect(1,0,1));
131  Box tbxyz = mfi.tilebox(IntVect(0,1,1));
132 
133  // We need a halo cell for terrain
134  bxcc.grow(IntVect(1,1,0));
135  tbxxy.grow(IntVect(1,1,0));
136  tbxxz.grow(IntVect(1,1,0));
137  tbxyz.grow(IntVect(1,1,0));
138 
139  if (bxcc.smallEnd(2) != domain.smallEnd(2)) {
140  bxcc.growLo(2,1);
141  tbxxy.growLo(2,1);
142  tbxxz.growLo(2,1);
143  tbxyz.growLo(2,1);
144  }
145 
146  if (bxcc.bigEnd(2) != domain.bigEnd(2)) {
147  bxcc.growHi(2,1);
148  tbxxy.growHi(2,1);
149  tbxxz.growHi(2,1);
150  tbxyz.growHi(2,1);
151  }
152 
153  // Expansion rate
154  Array4<Real> er_arr = expr->array(mfi);
155 
156  // Temporary storage for tiling/OMP
157  FArrayBox S11,S22,S33;
158  FArrayBox S12,S13,S23;
159  S11.resize( bxcc,1,The_Async_Arena()); S22.resize(bxcc,1,The_Async_Arena()); S33.resize(bxcc,1,The_Async_Arena());
160  S12.resize(tbxxy,1,The_Async_Arena()); S13.resize(tbxxz,1,The_Async_Arena()); S23.resize(tbxyz,1,The_Async_Arena());
161  Array4<Real> s11 = S11.array(); Array4<Real> s22 = S22.array(); Array4<Real> s33 = S33.array();
162  Array4<Real> s12 = S12.array(); Array4<Real> s13 = S13.array(); Array4<Real> s23 = S23.array();
163 
164  // Symmetric strain/stresses
165  Array4<Real> tau11 = Tau_lev[TauType::tau11]->array(mfi); Array4<Real> tau22 = Tau_lev[TauType::tau22]->array(mfi);
166  Array4<Real> tau33 = Tau_lev[TauType::tau33]->array(mfi); Array4<Real> tau12 = Tau_lev[TauType::tau12]->array(mfi);
167  Array4<Real> tau13 = Tau_lev[TauType::tau13]->array(mfi); Array4<Real> tau23 = Tau_lev[TauType::tau23]->array(mfi);
168 
169  // Strain magnitude
170  Array4<Real> SmnSmn_a;
171 
172  if (solverChoice.mesh_type == MeshType::StretchedDz) {
173  // Terrain non-symmetric terms
174  FArrayBox S21,S31,S32;
175  S21.resize(tbxxy,1,The_Async_Arena()); S31.resize(tbxxz,1,The_Async_Arena()); S32.resize(tbxyz,1,The_Async_Arena());
176  Array4<Real> s21 = S21.array(); Array4<Real> s31 = S31.array(); Array4<Real> s32 = S32.array();
177  Array4<Real> tau21 = Tau_lev[TauType::tau21]->array(mfi);
178  Array4<Real> tau31 = Tau_lev[TauType::tau31]->array(mfi);
179  Array4<Real> tau32 = Tau_lev[TauType::tau32]->array(mfi);
180 
181 
182  // *****************************************************************************
183  // Expansion rate compute terrain
184  // *****************************************************************************
185  {
186  BL_PROFILE("slow_rhs_making_er_S");
187  ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
188  {
189  Real mfsq = mf_mx(i,j,0)*mf_my(i,j,0);
190  er_arr(i,j,k) = (u(i+1, j , k )/mf_uy(i+1,j,0) - u(i, j, k)/mf_uy(i,j,0))*dxInv[0]*mfsq +
191  (v(i , j+1, k )/mf_vx(i,j+1,0) - v(i, j, k)/mf_vx(i,j,0))*dxInv[1]*mfsq +
192  (w(i , j , k+1) - w(i, j, k))/dz_ptr[k];
193  });
194  } // end profile
195 
196  // *****************************************************************************
197  // Strain tensor compute terrain
198  // *****************************************************************************
199  {
200  BL_PROFILE("slow_rhs_making_strain_S");
201  ComputeStrain_S(bxcc, tbxxy, tbxxz, tbxyz, domain,
202  u, v, w,
203  s11, s22, s33,
204  s12, s21,
205  s13, s31,
206  s23, s32,
207  stretched_dz_d, dxInv,
208  mf_mx, mf_ux, mf_vx,
209  mf_my, mf_uy, mf_vy, bc_ptr_h);
210  } // profile
211 
212  // Populate SmnSmn if using Deardorff or k-eqn RANS (used as diff src in post)
213  // and in the first RK stage (TKE tendencies constant for nrk>0, following WRF)
214  if ((nrk==0) && (need_SmnSmn)) {
215  SmnSmn_a = SmnSmn->array(mfi);
216  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
217  {
218  SmnSmn_a(i,j,k) = ComputeSmnSmn(i,j,k,s11,s22,s33,s12,s13,s23);
219  });
220  }
221 
222  // *****************************************************************************
223  // Stress tensor compute terrain
224  // *****************************************************************************
225  {
226  BL_PROFILE("slow_rhs_making_stress_T");
227 
228  // Remove Halo cells just for tau_ij comps
229  tbxxy.grow(IntVect(-1,-1,0));
230  tbxxz.grow(IntVect(-1,-1,0));
231  tbxyz.grow(IntVect(-1,-1,0));
232 
233  if (!l_use_turb) {
234  ComputeStressConsVisc_S(bxcc, tbxxy, tbxxz, tbxyz, mu_eff,
235  cell_data,
236  s11, s22, s33,
237  s12, s21,
238  s13, s31,
239  s23, s32,
240  er_arr, stretched_dz_d, dxInv,
241  mf_mx, mf_ux, mf_vx,
242  mf_my, mf_uy, mf_vy);
243  } else {
244  ComputeStressVarVisc_S(bxcc, tbxxy, tbxxz, tbxyz, mu_eff, mu_turb,
245  cell_data,
246  s11, s22, s33,
247  s12, s21,
248  s13, s31,
249  s23, s32,
250  er_arr, stretched_dz_d, dxInv,
251  mf_mx, mf_ux, mf_vx,
252  mf_my, mf_uy, mf_vy);
253  }
254 
255  // Remove halo cells from tau_ii but extend across valid_box bdry
256  bxcc.grow(IntVect(-1,-1,0));
257  if (bxcc.smallEnd(0) == valid_bx.smallEnd(0)) bxcc.growLo(0, 1);
258  if (bxcc.bigEnd(0) == valid_bx.bigEnd(0)) bxcc.growHi(0, 1);
259  if (bxcc.smallEnd(1) == valid_bx.smallEnd(1)) bxcc.growLo(1, 1);
260  if (bxcc.bigEnd(1) == valid_bx.bigEnd(1)) bxcc.growHi(1, 1);
261 
262  // Copy from temp FABs back to tau
263  ParallelFor(bxcc,
264  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
265  tau11(i,j,k) = s11(i,j,k);
266  tau22(i,j,k) = s22(i,j,k);
267  tau33(i,j,k) = s33(i,j,k);
268  });
269 
270  ParallelFor(tbxxy, tbxxz, tbxyz,
271  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
272  tau12(i,j,k) = s12(i,j,k);
273  tau21(i,j,k) = s21(i,j,k);
274  },
275  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
276  tau13(i,j,k) = s13(i,j,k);
277  tau31(i,j,k) = s31(i,j,k);
278  },
279  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
280  tau23(i,j,k) = s23(i,j,k);
281  tau32(i,j,k) = s32(i,j,k);
282  });
283  } // end profile
284 
285  } else if (l_use_terrain_fitted_coords) {
286 
287  // Terrain non-symmetric terms
288  FArrayBox S21,S31,S32;
289  S21.resize(tbxxy,1,The_Async_Arena()); S31.resize(tbxxz,1,The_Async_Arena()); S32.resize(tbxyz,1,The_Async_Arena());
290  Array4<Real> s21 = S21.array(); Array4<Real> s31 = S31.array(); Array4<Real> s32 = S32.array();
291  Array4<Real> tau21 = Tau_lev[TauType::tau21]->array(mfi);
292  Array4<Real> tau31 = Tau_lev[TauType::tau31]->array(mfi);
293  Array4<Real> tau32 = Tau_lev[TauType::tau32]->array(mfi);
294 
295 
296  // *****************************************************************************
297  // Expansion rate compute terrain
298  // *****************************************************************************
299  {
300  BL_PROFILE("slow_rhs_making_er_T");
301  Box gbxo = surroundingNodes(bxcc,2);
302 
303  // We make a temporary container for contravariant velocity Omega here
304  // -- it is only used to compute er_arr below
305  FArrayBox Omega;
306  Omega.resize(gbxo,1,The_Async_Arena());
307 
308  // First create Omega using velocity (not momentum)
309  Array4<Real> omega_arr = Omega.array();
310  ParallelFor(gbxo, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
311  {
312  omega_arr(i,j,k) = (k == 0) ? 0. : OmegaFromW(i,j,k,w(i,j,k),u,v,
313  mf_ux,mf_vy,z_nd,dxInv);
314  });
315 
316  ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
317  {
318 
319  Real met_u_h_zeta_hi = Compute_h_zeta_AtIface(i+1, j , k, dxInv, z_nd);
320  Real met_u_h_zeta_lo = Compute_h_zeta_AtIface(i , j , k, dxInv, z_nd);
321 
322  Real met_v_h_zeta_hi = Compute_h_zeta_AtJface(i , j+1, k, dxInv, z_nd);
323  Real met_v_h_zeta_lo = Compute_h_zeta_AtJface(i , j , k, dxInv, z_nd);
324 
325  Real Omega_hi = omega_arr(i,j,k+1);
326  Real Omega_lo = omega_arr(i,j,k );
327 
328  Real mfsq = mf_mx(i,j,0)*mf_my(i,j,0);
329 
330  Real expansionRate = (u(i+1,j ,k)/mf_uy(i+1,j,0)*met_u_h_zeta_hi - u(i,j,k)/mf_uy(i,j,0)*met_u_h_zeta_lo)*dxInv[0]*mfsq +
331  (v(i ,j+1,k)/mf_vx(i,j+1,0)*met_v_h_zeta_hi - v(i,j,k)/mf_vx(i,j,0)*met_v_h_zeta_lo)*dxInv[1]*mfsq +
332  (Omega_hi - Omega_lo)*dxInv[2];
333 
334  er_arr(i,j,k) = expansionRate / detJ_arr(i,j,k);
335  });
336  } // end profile
337 
338  // *****************************************************************************
339  // Strain tensor compute terrain
340  // *****************************************************************************
341  {
342  BL_PROFILE("slow_rhs_making_strain_T");
343  ComputeStrain_T(bxcc, tbxxy, tbxxz, tbxyz, domain,
344  u, v, w,
345  s11, s22, s33,
346  s12, s21,
347  s13, s31,
348  s23, s32,
349  z_nd, detJ_arr, dxInv,
350  mf_mx, mf_ux, mf_vx,
351  mf_my, mf_uy, mf_vy, bc_ptr_h);
352  } // profile
353 
354  // Populate SmnSmn if using Deardorff or k-eqn RANS (used as diff src in post)
355  // and in the first RK stage (TKE tendencies constant for nrk>0, following WRF)
356  if ((nrk==0) && (need_SmnSmn)) {
357  SmnSmn_a = SmnSmn->array(mfi);
358  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
359  {
360  SmnSmn_a(i,j,k) = ComputeSmnSmn(i,j,k,
361  s11,s22,s33,
362  s12,s13,s23);
363  });
364  }
365 
366  // *****************************************************************************
367  // Stress tensor compute terrain
368  // *****************************************************************************
369  {
370  BL_PROFILE("slow_rhs_making_stress_T");
371 
372  // Remove Halo cells just for tau_ij comps
373  tbxxy.grow(IntVect(-1,-1,0));
374  tbxxz.grow(IntVect(-1,-1,0));
375  tbxyz.grow(IntVect(-1,-1,0));
376 
377  if (!l_use_turb) {
378  ComputeStressConsVisc_T(bxcc, tbxxy, tbxxz, tbxyz, mu_eff,
379  cell_data,
380  s11, s22, s33,
381  s12, s21,
382  s13, s31,
383  s23, s32,
384  er_arr, z_nd, detJ_arr, dxInv,
385  mf_mx, mf_ux, mf_vx,
386  mf_my, mf_uy, mf_vy);
387  } else {
388  ComputeStressVarVisc_T(bxcc, tbxxy, tbxxz, tbxyz, mu_eff, mu_turb,
389  cell_data,
390  s11, s22, s33,
391  s12, s21,
392  s13, s31,
393  s23, s32,
394  er_arr, z_nd, detJ_arr, dxInv,
395  mf_mx, mf_ux, mf_vx,
396  mf_my, mf_uy, mf_vy);
397  }
398 
399  // Remove halo cells from tau_ii but extend across valid_box bdry
400  bxcc.grow(IntVect(-1,-1,0));
401  if (bxcc.smallEnd(0) == valid_bx.smallEnd(0)) bxcc.growLo(0, 1);
402  if (bxcc.bigEnd(0) == valid_bx.bigEnd(0)) bxcc.growHi(0, 1);
403  if (bxcc.smallEnd(1) == valid_bx.smallEnd(1)) bxcc.growLo(1, 1);
404  if (bxcc.bigEnd(1) == valid_bx.bigEnd(1)) bxcc.growHi(1, 1);
405 
406  // Copy from temp FABs back to tau
407  ParallelFor(bxcc,
408  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
409  tau11(i,j,k) = s11(i,j,k);
410  tau22(i,j,k) = s22(i,j,k);
411  tau33(i,j,k) = s33(i,j,k);
412  });
413 
414  ParallelFor(tbxxy, tbxxz, tbxyz,
415  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
416  tau12(i,j,k) = s12(i,j,k);
417  tau21(i,j,k) = s21(i,j,k);
418  },
419  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
420  tau13(i,j,k) = s13(i,j,k);
421  tau31(i,j,k) = s31(i,j,k);
422  },
423  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
424  tau23(i,j,k) = s23(i,j,k);
425  tau32(i,j,k) = s32(i,j,k);
426  });
427  } // end profile
428 
429  } else {
430 
431  // *****************************************************************************
432  // Expansion rate compute no terrain
433  // *****************************************************************************
434  {
435  BL_PROFILE("slow_rhs_making_er_N");
436  ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
437  Real mfsq = mf_mx(i,j,0)*mf_my(i,j,0);
438  er_arr(i,j,k) = (u(i+1, j , k )/mf_uy(i+1,j,0) - u(i, j, k)/mf_uy(i,j,0))*dxInv[0]*mfsq +
439  (v(i , j+1, k )/mf_vx(i,j+1,0) - v(i, j, k)/mf_vx(i,j,0))*dxInv[1]*mfsq +
440  (w(i , j , k+1) - w(i, j, k))*dxInv[2];
441  });
442  } // end profile
443 
444 
445  // *****************************************************************************
446  // Strain tensor compute no terrain
447  // *****************************************************************************
448  {
449  BL_PROFILE("slow_rhs_making_strain_N");
450  ComputeStrain_N(bxcc, tbxxy, tbxxz, tbxyz, domain,
451  u, v, w,
452  s11, s22, s33,
453  s12, s13, s23,
454  dxInv,
455  mf_mx, mf_ux, mf_vx,
456  mf_my, mf_uy, mf_vy, bc_ptr_h);
457  } // end profile
458 
459  // Populate SmnSmn if using Deardorff or k-eqn RANS (used as diff src in post)
460  // and in the first RK stage (TKE tendencies constant for nrk>0, following WRF)
461  if ((nrk==0) && (need_SmnSmn)) {
462  SmnSmn_a = SmnSmn->array(mfi);
463  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
464  {
465  SmnSmn_a(i,j,k) = ComputeSmnSmn(i,j,k,
466  s11,s22,s33,
467  s12,s13,s23);
468  });
469  }
470 
471  // *****************************************************************************
472  // Stress tensor compute no terrain
473  // *****************************************************************************
474  {
475  BL_PROFILE("slow_rhs_making_stress_N");
476 
477  // Remove Halo cells just for tau_ij comps
478  tbxxy.grow(IntVect(-1,-1,0));
479  tbxxz.grow(IntVect(-1,-1,0));
480  tbxyz.grow(IntVect(-1,-1,0));
481  if (tbxxy.smallEnd(2) > domlo_z) {
482  tbxxy.growLo(2,-1);
483  tbxxz.growLo(2,-1);
484  tbxyz.growLo(2,-1);
485  }
486  if (tbxxy.bigEnd(2) < domhi_z) {
487  tbxxy.growHi(2,-1);
488  tbxxz.growHi(2,-1);
489  tbxyz.growHi(2,-1);
490  }
491 
492  if (!l_use_turb) {
493  ComputeStressConsVisc_N(bxcc, tbxxy, tbxxz, tbxyz, mu_eff,
494  cell_data,
495  s11, s22, s33,
496  s12, s13, s23,
497  er_arr);
498  } else {
499  ComputeStressVarVisc_N(bxcc, tbxxy, tbxxz, tbxyz, mu_eff, mu_turb,
500  cell_data,
501  s11, s22, s33,
502  s12, s13, s23,
503  er_arr);
504  }
505 
506  // Remove halo cells from tau_ii but extend across valid_box bdry
507  bxcc.grow(IntVect(-1,-1,0));
508  if (bxcc.smallEnd(0) == valid_bx.smallEnd(0)) bxcc.growLo(0, 1);
509  if (bxcc.bigEnd(0) == valid_bx.bigEnd(0)) bxcc.growHi(0, 1);
510  if (bxcc.smallEnd(1) == valid_bx.smallEnd(1)) bxcc.growLo(1, 1);
511  if (bxcc.bigEnd(1) == valid_bx.bigEnd(1)) bxcc.growHi(1, 1);
512 
513  // Copy from temp FABs back to tau
514  ParallelFor(bxcc,
515  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
516  tau11(i,j,k) = s11(i,j,k);
517  tau22(i,j,k) = s22(i,j,k);
518  tau33(i,j,k) = s33(i,j,k);
519  });
520  ParallelFor(tbxxy, tbxxz, tbxyz,
521  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
522  tau12(i,j,k) = s12(i,j,k);
523  },
524  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
525  tau13(i,j,k) = s13(i,j,k);
526  },
527  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
528  tau23(i,j,k) = s23(i,j,k);
529  });
530  } // end profile
531  } // l_use_terrain_fitted_coords
532  } // MFIter
533  } // l_use_diff
534 }
void ComputeStrain_N(Box bxcc, Box tbxxy, Box tbxxz, Box tbxyz, Box domain, const Array4< const Real > &u, const Array4< const Real > &v, const Array4< const Real > &w, Array4< Real > &tau11, Array4< Real > &tau22, Array4< Real > &tau33, Array4< Real > &tau12, Array4< Real > &tau13, Array4< Real > &tau23, 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, const BCRec *bc_ptr)
Definition: ERF_ComputeStrain_N.cpp:28
void ComputeStrain_S(Box bxcc, Box tbxxy, Box tbxxz, Box tbxyz, Box domain, const Array4< const Real > &u, const Array4< const Real > &v, const Array4< const Real > &w, 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 Gpu::DeviceVector< Real > &stretched_dz_d, 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, const BCRec *bc_ptr)
Definition: ERF_ComputeStrain_S.cpp:36
void ComputeStrain_T(Box bxcc, Box tbxxy, Box tbxxz, Box tbxyz, Box domain, const Array4< const Real > &u, const Array4< const Real > &v, const Array4< const Real > &w, 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 > &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, const BCRec *bc_ptr)
Definition: ERF_ComputeStrain_T.cpp:36
void ComputeStressVarVisc_N(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 > &tau23, const Array4< const Real > &er_arr)
Definition: ERF_ComputeStress_N.cpp:105
void ComputeStressConsVisc_N(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 > &tau23, const Array4< const Real > &er_arr)
Definition: ERF_ComputeStress_N.cpp:23
void ComputeStressVarVisc_S(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 Gpu::DeviceVector< Real > &stretched_dz_d, 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)
Definition: ERF_ComputeStress_S.cpp:150
void ComputeStressConsVisc_S(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 Gpu::DeviceVector< Real > &stretched_dz_d, 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)
Definition: ERF_ComputeStress_S.cpp:29
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)
Definition: ERF_ComputeStress_T.cpp:29
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)
Definition: ERF_ComputeStress_T.cpp:336
@ 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
@ v_x
Definition: ERF_DataStruct.H:22
@ u_y
Definition: ERF_DataStruct.H:23
@ v_y
Definition: ERF_DataStruct.H:23
@ m_y
Definition: ERF_DataStruct.H:23
@ u_x
Definition: ERF_DataStruct.H:22
@ m_x
Definition: ERF_DataStruct.H:22
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real ComputeSmnSmn(int &i, int &j, int &k, const amrex::Array4< amrex::Real const > &tau11, const amrex::Array4< amrex::Real const > &tau22, const amrex::Array4< amrex::Real const > &tau33, const amrex::Array4< amrex::Real const > &tau12, const amrex::Array4< amrex::Real const > &tau13, const amrex::Array4< amrex::Real const > &tau23)
Definition: ERF_EddyViscosity.H:40
amrex::Real Real
Definition: ERF_ShocInterface.H:19
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real OmegaFromW(int &i, int &j, int &k, amrex::Real w, const amrex::Array4< const amrex::Real > &u_arr, const amrex::Array4< const amrex::Real > &v_arr, const amrex::Array4< const amrex::Real > &mf_u, const amrex::Array4< const amrex::Real > &mf_v, const amrex::Array4< const amrex::Real > &z_nd, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &dxInv)
Definition: ERF_TerrainMetrics.H:415
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_zeta_AtIface(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:96
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_zeta_AtJface(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:139
AMREX_FORCE_INLINE amrex::IntVect TileNoZ()
Definition: ERF_TileNoZ.H:11
@ cons
Definition: ERF_IndexDefines.H:158
@ xvel
Definition: ERF_IndexDefines.H:141
@ zvel
Definition: ERF_IndexDefines.H:143
@ yvel
Definition: ERF_IndexDefines.H:142
Definition: ERF_DiffStruct.H:19
amrex::Real rho0_trans
Definition: ERF_DiffStruct.H:91
MolecDiffType molec_diff_type
Definition: ERF_DiffStruct.H:84
amrex::Real dynamic_viscosity
Definition: ERF_DiffStruct.H:96
static MeshType mesh_type
Definition: ERF_DataStruct.H:815
DiffChoice diffChoice
Definition: ERF_DataStruct.H:824
amrex::Vector< TurbChoice > turbChoice
Definition: ERF_DataStruct.H:826
static TerrainType terrain_type
Definition: ERF_DataStruct.H:806
Definition: ERF_TurbStruct.H:41
PBLType pbl_type
Definition: ERF_TurbStruct.H:375
RANSType rans_type
Definition: ERF_TurbStruct.H:372
LESType les_type
Definition: ERF_TurbStruct.H:330
bool use_kturb
Definition: ERF_TurbStruct.H:381

Referenced by erf_slow_rhs_pre().

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