ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
Functions
ERF_DiffusionSrcForMom_EB.cpp File Reference
#include <AMReX.H>
#include <AMReX_EB_Slopes_K.H>
#include <ERF_EB.H>
#include <ERF_Diffusion.H>
#include <ERF_IndexDefines.H>
#include <ERF_EBSlopes.H>
#include <ERF_DiffStruct.H>
#include <ERF_EBStruct.H>
Include dependency graph for ERF_DiffusionSrcForMom_EB.cpp:

Functions

AMREX_GPU_DEVICE AMREX_FORCE_INLINE void compute_tangent_vectors (Real nx, Real ny, Real nz, Real &tbx_x, Real &tbx_y, Real &tbx_z, Real &tby_x, Real &tby_y, Real &tby_z)
 
void DiffusionSrcForMom_EB (const MFIter &mfi, [[maybe_unused]] const Box &domain, const Box &bxx, const Box &bxy, const Box &bxz, const Array4< Real > &rho_u_rhs, const Array4< Real > &rho_v_rhs, const Array4< Real > &rho_w_rhs, const Array4< const Real > &u_arr, const Array4< const Real > &v_arr, const Array4< const Real > &w_arr, const Array4< const Real > &tau11, const Array4< const Real > &tau22, const Array4< const Real > &tau33, const Array4< const Real > &tau12, const Array4< const Real > &tau13, const Array4< const Real > &tau23, const Array4< const Real > &u_tau_eb13, const Array4< const Real > &u_tau_eb23, const Array4< const Real > &v_tau_eb13, const Array4< const Real > &v_tau_eb23, const Array4< const Real > &w_tau_eb13, const Array4< const Real > &w_tau_eb23, const Real *dx_arr, 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 SolverChoice &solverChoice, const eb_ &ebfact, [[maybe_unused]] const BCRec *d_bcrec_ptr)
 

Function Documentation

◆ compute_tangent_vectors()

AMREX_GPU_DEVICE AMREX_FORCE_INLINE void compute_tangent_vectors ( Real  nx,
Real  ny,
Real  nz,
Real &  tbx_x,
Real &  tbx_y,
Real &  tbx_z,
Real &  tby_x,
Real &  tby_y,
Real &  tby_z 
)

Compute the tangent directions at an EB boundary given the unit normal vector.

t_bx and t_by are the unit vectors along the projections of e_x and e_y onto the tangent plane. They are the directions in which the stored EB surface stresses tau_eb13 and tau_eb23 act, since those are built from the x- and y-components of the tangential velocity (see ERF_EBMOSTStress.H). For a horizontal surface they are exactly e_x and e_y. Note that they are individually unit but are not orthogonal to each other unless nx*ny is zero – t_bx . t_by = -nx*ny.

e_x has no projection onto the tangent plane when n is parallel to e_x, and likewise for e_y, so one of the two is undefined for an axis-aligned normal such as (1,0,0) (issue #3533). That limit is genuinely two-sided – approaching n = (1,0,0) from nz > 0 and from nz < 0 gives t_bx of -e_z and +e_z – so there is no continuous choice, and the degenerate vector is returned as zero. That keeps the momentum RHS finite, and it keeps the component that does have a limit continuous: the x-component of t_bx is sqrt(1-nx*nx), which tends to zero here as well. It also avoids driving w-momentum with tau_eb13, an x-associated stress, along an arbitrarily signed direction. Only nx and ny can be degenerate, and never both at once.

Parameters
[in]nxx-component of the normal vector
[in]nyy-component of the normal vector
[in]nzz-component of the normal vector
[out]tbx_xx-component of first tangent vector
[out]tbx_yy-component of first tangent vector
[out]tbx_zz-component of first tangent vector
[out]tby_xx-component of second tangent vector
[out]tby_yy-component of second tangent vector
[out]tby_zz-component of second tangent vector
46 {
47  // Below this the projection is all round-off and its direction is meaningless.
48  // The squared norms are formed as ny^2+nz^2 and nx^2+nz^2 rather than as the
49  // algebraically equal 1-nx^2 and 1-ny^2 to avoid cancellation for a near-axis normal
50  // The squared norm of a unit normal carries an absolute roundoff of order
51  // real_eps, so the floor is a (generous) multiple of it: 2.2e-12 in double,
52  // 1.2e-03 in single, i.e. a normal within ~2 degrees of the axis.
53  constexpr Real tol = Real(1.e4)*real_eps;
54 
55  // x-tangential vector: t_bx = (e_x - (e_x · n)n) / ||e_x - (e_x · n)n||
56  // e_x = (1,0,0), so e_x · n = nx
57  Real tbx_norm2 = ny*ny + nz*nz;
58  if (tbx_norm2 > tol) {
59  Real tbx_norm_inv = one / std::sqrt(tbx_norm2);
60  tbx_x = ( one - nx * nx) * tbx_norm_inv;
61  tbx_y = ( - nx * ny) * tbx_norm_inv;
62  tbx_z = ( - nx * nz) * tbx_norm_inv;
63  } else {
64  tbx_x = zero;
65  tbx_y = zero;
66  tbx_z = zero;
67  }
68 
69  // y-tangential vector: t_by = (e_y - (e_y · n)n) / ||e_y - (e_y · n)n||
70  // e_y = (0,1,0), so e_y · n = ny
71  Real tby_norm2 = nx*nx + nz*nz;
72  if (tby_norm2 > tol) {
73  Real tby_norm_inv = one / std::sqrt(tby_norm2);
74  tby_x = ( - ny * nx) * tby_norm_inv;
75  tby_y = ( one - ny * ny) * tby_norm_inv;
76  tby_z = ( - ny * nz) * tby_norm_inv;
77  } else {
78  tby_x = zero;
79  tby_y = zero;
80  tby_z = zero;
81  }
82 }
nx
const int nx
Definition: ERF_InitCustomPertVels_CloudChamber.H:14
ny
const int ny
Definition: ERF_InitCustomPertVels_CloudChamber.H:15
real_eps
constexpr amrex::Real real_eps
Definition: ERF_NumericalConstants.H:27
one
constexpr amrex::Real one
Definition: ERF_NumericalConstants.H:30
zero
constexpr amrex::Real zero
Definition: ERF_NumericalConstants.H:29
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:19

Referenced by DiffusionSrcForMom_EB().

Here is the caller graph for this function:

◆ DiffusionSrcForMom_EB()

void DiffusionSrcForMom_EB ( const MFIter &  mfi,
[[maybe_unused] ] const Box &  domain,
const Box &  bxx,
const Box &  bxy,
const Box &  bxz,
const Array4< Real > &  rho_u_rhs,
const Array4< Real > &  rho_v_rhs,
const Array4< Real > &  rho_w_rhs,
const Array4< const Real > &  u_arr,
const Array4< const Real > &  v_arr,
const Array4< const Real > &  w_arr,
const Array4< const Real > &  tau11,
const Array4< const Real > &  tau22,
const Array4< const Real > &  tau33,
const Array4< const Real > &  tau12,
const Array4< const Real > &  tau13,
const Array4< const Real > &  tau23,
const Array4< const Real > &  u_tau_eb13,
const Array4< const Real > &  u_tau_eb23,
const Array4< const Real > &  v_tau_eb13,
const Array4< const Real > &  v_tau_eb23,
const Array4< const Real > &  w_tau_eb13,
const Array4< const Real > &  w_tau_eb23,
const Real *  dx_arr,
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 SolverChoice &  solverChoice,
const eb_ &  ebfact,
[[maybe_unused] ] const BCRec *  d_bcrec_ptr 
)

Function for computing the momentum RHS for diffusion operator without terrain.

Parameters
[in]mfiMultiFab Iterator
[in]domaincomputational domain
[in]bxxnodal x box for x-mom
[in]bxynodal y box for y-mom
[in]bxznodal z box for z-mom
[out]rho_u_rhsRHS for x-mom
[out]rho_v_rhsRHS for y-mom
[out]rho_w_rhsRHS for z-mom
[in]u_arrx-direction velocity
[in]v_arry-direction velocity
[in]w_arrz-direction velocity
[in]tau1111 stress
[in]tau2222 stress
[in]tau3333 stress
[in]tau1212 stress
[in]tau1313 stress
[in]tau2323 stress
[in]u_tau_eb13EB tangential stress for x-momentum on z-faces
[in]u_tau_eb23EB tangential stress for x-momentum on y-faces
[in]v_tau_eb13EB tangential stress for y-momentum on z-faces
[in]v_tau_eb23EB tangential stress for y-momentum on x-faces
[in]w_tau_eb13EB tangential stress for z-momentum on y-faces
[in]w_tau_eb23EB tangential stress for z-momentum on x-faces
[in]dx_arrcell size array
[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]solverChoicecontainer with diffusion parameters
[in]ebfactEB factories for cell- and face-centered variables
[in]d_bcrec_ptrboundary condition records on device
155 {
156  BL_PROFILE_VAR("DiffusionSrcForMom_EB()",DiffusionSrcForMom_EB);
157 
158  DiffChoice dc = solverChoice.diffChoice;
159  const bool l_use_constAlpha = ( dc.molec_diff_type == MolecDiffType::ConstantAlpha );
160  Real mu_eff = (l_use_constAlpha) ? two * dc.dynamic_viscosity / dc.rho0_trans
161  : two * dc.dynamic_viscosity;
162 
163  auto dxinv = dxInv[0], dyinv = dxInv[1], dzinv = dxInv[2];
164  Real dx = dx_arr[0], dy = dx_arr[1], dz = dx_arr[2];
165  Real vol = dx * dy * dz;
166 
167  EBChoice ebChoice = solverChoice.ebChoice;
168  const bool l_no_slip = (ebChoice.eb_boundary_type == EBBoundaryType::NoSlipWall);
169  const bool l_surface_layer = (ebChoice.eb_boundary_type == EBBoundaryType::SurfaceLayer);
170  const bool l_constraint_x = solverChoice.diffChoice.eb_diff_constraint_x;
171  const bool l_constraint_y = solverChoice.diffChoice.eb_diff_constraint_y;
172  const bool l_constraint_z = solverChoice.diffChoice.eb_diff_constraint_z;
173 
174  // EB u-factory
175  const auto* u_factory = ebfact.get_u_const_factory();
176  Array4<const EBCellFlag> u_cellflg = u_factory->getMultiEBCellFlagFab()[mfi].const_array();
177  Array4<const Real > u_volfrac = u_factory->getVolFrac().const_array(mfi);
178  Array4<const Real > u_volcent{};
179  Array4<const Real > u_afrac_x{};
180  Array4<const Real > u_afrac_y{};
181  Array4<const Real > u_afrac_z{};
182  Array4<const Real > u_bcent{};
183  Array4<const Real > u_bnorm{};
184  FabType u_type = u_factory->getMultiEBCellFlagFab()[mfi].getType();
185  if (u_type == FabType::singlevalued) {
186  u_volcent = u_factory->getCentroid().const_array(mfi);
187  u_afrac_x = u_factory->getAreaFrac()[0]->const_array(mfi);
188  u_afrac_y = u_factory->getAreaFrac()[1]->const_array(mfi);
189  u_afrac_z = u_factory->getAreaFrac()[2]->const_array(mfi);
190  u_bcent = u_factory->getBndryCent().const_array(mfi);
191  u_bnorm = u_factory->getBndryNormal().const_array(mfi);
192  }
193 
194 
195  // EB v-factory
196  const auto* v_factory = ebfact.get_v_const_factory();
197  Array4<const EBCellFlag> v_cellflg = v_factory->getMultiEBCellFlagFab()[mfi].const_array();
198  Array4<const Real > v_volfrac = v_factory->getVolFrac().const_array(mfi);
199  Array4<const Real > v_volcent{};
200  Array4<const Real > v_afrac_x{};
201  Array4<const Real > v_afrac_y{};
202  Array4<const Real > v_afrac_z{};
203  Array4<const Real > v_bcent{};
204  Array4<const Real > v_bnorm{};
205  FabType v_type = v_factory->getMultiEBCellFlagFab()[mfi].getType();
206  if (v_type == FabType::singlevalued) {
207  v_volcent = v_factory->getCentroid().const_array(mfi);
208  v_afrac_x = v_factory->getAreaFrac()[0]->const_array(mfi);
209  v_afrac_y = v_factory->getAreaFrac()[1]->const_array(mfi);
210  v_afrac_z = v_factory->getAreaFrac()[2]->const_array(mfi);
211  v_bcent = v_factory->getBndryCent().const_array(mfi);
212  v_bnorm = v_factory->getBndryNormal().const_array(mfi);
213  }
214 
215  // EB w-factory
216  const auto* w_factory = ebfact.get_w_const_factory();
217  Array4<const EBCellFlag> w_cellflg = w_factory->getMultiEBCellFlagFab()[mfi].const_array();
218  Array4<const Real > w_volfrac = w_factory->getVolFrac().const_array(mfi);
219  Array4<const Real > w_volcent{};
220  Array4<const Real > w_afrac_x{};
221  Array4<const Real > w_afrac_y{};
222  Array4<const Real > w_afrac_z{};
223  Array4<const Real > w_bcent{};
224  Array4<const Real > w_bnorm{};
225  FabType w_type = w_factory->getMultiEBCellFlagFab()[mfi].getType();
226  if (w_type == FabType::singlevalued) {
227  w_volcent = w_factory->getCentroid().const_array(mfi);
228  w_afrac_x = w_factory->getAreaFrac()[0]->const_array(mfi);
229  w_afrac_y = w_factory->getAreaFrac()[1]->const_array(mfi);
230  w_afrac_z = w_factory->getAreaFrac()[2]->const_array(mfi);
231  w_bcent = w_factory->getBndryCent().const_array(mfi);
232  w_bnorm = w_factory->getBndryNormal().const_array(mfi);
233  }
234 
235  // x-momentum
236  if (u_type == FabType::regular) {
237 
238  ParallelFor(bxx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
239 
240  Real mfsq = mf_ux(i,j,0) * mf_uy(i,j,0);
241 
242  Real diffContrib = ( (tau11(i , j , k ) - tau11(i-1, j , k ) ) * dxinv * mfsq
243  + (tau12(i , j+1, k ) - tau12(i , j , k ) ) * dyinv * mfsq
244  + (tau13(i , j , k+1) - tau13(i , j , k ) ) * dzinv );
245  diffContrib /= u_volfrac(i,j,k);
246 
247  rho_u_rhs(i,j,k) -= diffContrib;
248  });
249 
250  } else if (u_type == FabType::singlevalued) {
251 
252  ParallelFor(bxx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
253 
254  if (u_volfrac(i,j,k)>zero) {
255 
256  // Inv Jacobian
257  Real mfsq = mf_ux(i,j,0) * mf_uy(i,j,0);
258 
259  Real diffContrib = ( (tau11(i , j , k ) * u_afrac_x(i+1,j ,k )
260  - tau11(i-1, j , k ) * u_afrac_x(i ,j ,k ) ) * dxinv * mfsq
261  + (tau12(i , j+1, k ) * u_afrac_y(i ,j+1,k )
262  - tau12(i , j , k ) * u_afrac_y(i ,j ,k ) ) * dyinv * mfsq
263  + (tau13(i , j , k+1) * u_afrac_z(i ,j ,k+1)
264  - tau13(i , j , k ) * u_afrac_z(i ,j ,k )) * dzinv );
265  diffContrib /= u_volfrac(i,j,k);
266 
267  rho_u_rhs(i,j,k) -= diffContrib;
268 
269  if (!l_constraint_x && u_cellflg(i,j,k).isSingleValued()) {
270 
271  Real axm = u_afrac_x(i ,j ,k );
272  Real axp = u_afrac_x(i+1,j ,k );
273  Real aym = u_afrac_y(i ,j ,k );
274  Real ayp = u_afrac_y(i ,j+1,k );
275  Real azm = u_afrac_z(i ,j ,k );
276  Real azp = u_afrac_z(i ,j ,k+1);
277 
278  Real adx = (axm-axp) * dy * dz;
279  Real ady = (aym-ayp) * dx * dz;
280  Real adz = (azm-azp) * dx * dy;
281 
282  Real barea = std::sqrt(adx*adx + ady*ady + adz*adz);
283 
284  Real dudn = zero;
285 
286  if (l_no_slip || l_surface_layer) {
287 
288  RealVect bcent_eb {u_bcent(i,j,k,0), u_bcent(i,j,k,1), u_bcent(i,j,k,2)};
289 
290  Real Dirichlet_u {zero};
291  Real Dirichlet_v {zero};
292  Real Dirichlet_w {zero};
293 
294  Real nx = u_bnorm(i,j,k,0);
295  Real ny = u_bnorm(i,j,k,1);
296  Real nz = u_bnorm(i,j,k,2);
297 
298  if (l_surface_layer) {
299 
300  // Average v and w onto the x-face
301  Real velx = u_arr(i,j,k);
302  Real vely = (v_volfrac(i-1,j ,k) * v_arr(i-1,j ,k) + v_volfrac(i,j ,k) * v_arr(i,j ,k)
303  + v_volfrac(i-1,j+1,k) * v_arr(i-1,j+1,k) + v_volfrac(i,j+1,k) * v_arr(i,j+1,k))
304  / (v_volfrac(i-1,j,k) + v_volfrac(i,j,k) + v_volfrac(i-1,j+1,k) + v_volfrac(i,j+1,k));
305 
306  Real velz = (w_volfrac(i-1,j,k ) * w_arr(i-1,j,k ) + w_volfrac(i,j,k ) * w_arr(i,j,k )
307  + w_volfrac(i-1,j,k+1) * w_arr(i-1,j,k+1) + w_volfrac(i,j,k+1) * w_arr(i,j,k+1))
308  / (w_volfrac(i-1,j,k) + w_volfrac(i,j,k) + w_volfrac(i-1,j,k+1) + w_volfrac(i,j,k+1));
309 
310  // Impose tangential velocity as Dirichlet condition
311  Real v_dot_n = velx * nx + vely * ny + velz * nz;
312  Dirichlet_u = velx - v_dot_n * nx;
313  Dirichlet_v = vely - v_dot_n * ny;
314  Dirichlet_w = velz - v_dot_n * nz;
315  }
316 
317  GpuArray<Real,AMREX_SPACEDIM> slopes_u;
318  GpuArray<Real,AMREX_SPACEDIM> slopes_v;
319  GpuArray<Real,AMREX_SPACEDIM> slopes_w;
320 
321  slopes_u = erf_calc_slopes_eb_Dirichlet ( dx, dy, dz, i, j, k, bcent_eb, Dirichlet_u, u_arr, u_volcent, u_cellflg);
322  slopes_v = erf_calc_slopes_eb_Dirichlet_staggered( Vars::xvel, Vars::yvel, dx, dy, dz, i, j, k, bcent_eb, Dirichlet_v, v_arr, v_volcent, v_cellflg);
323  slopes_w = erf_calc_slopes_eb_Dirichlet_staggered( Vars::xvel, Vars::zvel, dx, dy, dz, i, j, k, bcent_eb, Dirichlet_w, w_arr, w_volcent, w_cellflg);
324 
325  Real dudx = slopes_u[0];
326  Real dudy = slopes_u[1];
327  Real dudz = slopes_u[2];
328  Real dvdx = slopes_v[0];
329  Real dvdy = slopes_v[1];
330  Real dvdz = slopes_v[2];
331  Real dwdx = slopes_w[0];
332  Real dwdy = slopes_w[1];
333  Real dwdz = slopes_w[2];
334 
335  Real tau11_eb = ( dudx - ( dudx + dvdy + dwdz ) / three );
336  Real tau12_eb = myhalf * (dudy + dvdx);
337  Real tau13_eb = myhalf * (dudz + dwdx);
338 
339  if (l_no_slip) {
340 
341  dudn = - mu_eff * (nx * tau11_eb + ny * tau12_eb + nz * tau13_eb);
342 
343  } else if (l_surface_layer) {
344 
345  Real tbx_x, tbx_y, tbx_z, tby_x, tby_y, tby_z;
346  compute_tangent_vectors(nx, ny, nz, tbx_x, tbx_y, tbx_z, tby_x, tby_y, tby_z);
347 
348  Real tau22_eb = ( dvdy - ( dudx + dvdy + dwdz ) / three );
349  Real tau33_eb = ( dwdz - ( dudx + dvdy + dwdz ) / three );
350  Real tau23_eb = myhalf * (dvdz + dwdy);
351 
352  Real tauzz = mu_eff * ( nx*nx*tau11_eb + ny*ny*tau22_eb + nz*nz*tau33_eb
353  + two * (nx*ny*tau12_eb + ny*nz*tau23_eb + nx*nz*tau13_eb ));
354 
355  dudn = - tbx_x * u_tau_eb13(i,j,k) - tby_x * u_tau_eb23(i,j,k) - nx * tauzz;
356  }
357  }
358 
359  rho_u_rhs(i,j,k) -= barea * dudn / (vol * u_volfrac(i,j,k));
360  }
361  }
362 
363  });
364  } // u_type
365 
366  // y-momentum
367  if (v_type == FabType::regular) {
368 
369  ParallelFor(bxy, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
370 
371  Real mfsq = mf_vx(i,j,0) * mf_vy(i,j,0);
372 
373  Real diffContrib = ( (tau12(i+1, j , k ) - tau12(i , j , k ) ) * dxinv * mfsq
374  + (tau22(i , j , k ) - tau22(i , j-1, k ) ) * dyinv * mfsq
375  + (tau23(i , j , k+1) - tau23(i , j , k ) ) * dzinv );
376  diffContrib /= v_volfrac(i,j,k);
377 
378  rho_v_rhs(i,j,k) -= diffContrib;
379  });
380  } else if (v_type == FabType::singlevalued) {
381 
382  ParallelFor(bxy, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
383 
384  if (v_volfrac(i,j,k)>zero) {
385 
386  // Inv Jacobian
387  Real mfsq = mf_vx(i,j,0) * mf_vy(i,j,0);
388 
389  Real diffContrib = ( (tau12(i+1, j , k ) * v_afrac_x(i+1,j ,k )
390  - tau12(i , j , k ) * v_afrac_x(i ,j ,k ) ) * dxinv * mfsq
391  + (tau22(i , j , k ) * v_afrac_y(i ,j+1,k )
392  - tau22(i , j-1, k ) * v_afrac_y(i ,j ,k ) ) * dyinv * mfsq
393  + (tau23(i , j , k+1) * v_afrac_z(i ,j ,k+1)
394  - tau23(i , j , k ) * v_afrac_z(i ,j ,k ) ) * dzinv );
395  diffContrib /= v_volfrac(i,j,k);
396 
397  rho_v_rhs(i,j,k) -= diffContrib;
398 
399  if (!l_constraint_y && v_cellflg(i,j,k).isSingleValued()) {
400 
401  Real axm = v_afrac_x(i ,j ,k );
402  Real axp = v_afrac_x(i+1,j ,k );
403  Real aym = v_afrac_y(i ,j ,k );
404  Real ayp = v_afrac_y(i ,j+1,k );
405  Real azm = v_afrac_z(i ,j ,k );
406  Real azp = v_afrac_z(i ,j ,k+1);
407 
408  Real adx = (axm-axp) * dy * dz;
409  Real ady = (aym-ayp) * dx * dz;
410  Real adz = (azm-azp) * dx * dy;
411 
412  Real barea = std::sqrt(adx*adx + ady*ady + adz*adz);
413 
414  Real dvdn = 0.0;
415 
416  if (l_no_slip || l_surface_layer) {
417 
418  RealVect bcent_eb {v_bcent(i,j,k,0), v_bcent(i,j,k,1), v_bcent(i,j,k,2)};
419 
420  Real Dirichlet_u {zero};
421  Real Dirichlet_v {zero};
422  Real Dirichlet_w {zero};
423 
424  Real nx = v_bnorm(i,j,k,0);
425  Real ny = v_bnorm(i,j,k,1);
426  Real nz = v_bnorm(i,j,k,2);
427 
428  if (l_surface_layer) {
429 
430  // Average u and w onto the y-face
431  Real velx = (u_volfrac(i ,j-1,k) * u_arr(i ,j-1,k) + u_volfrac(i+1,j-1,k) * u_arr(i+1,j-1,k)
432  + u_volfrac(i+1,j ,k) * u_arr(i+1,j ,k) + u_volfrac(i ,j ,k) * u_arr(i ,j ,k))
433  / (u_volfrac(i,j-1,k) + u_volfrac(i+1,j-1,k) + u_volfrac(i+1,j,k) + u_volfrac(i,j,k));
434  Real vely = v_arr(i,j,k);
435  Real velz = (w_volfrac(i,j-1,k ) * w_arr(i,j-1,k ) + w_volfrac(i,j,k ) * w_arr(i,j,k )
436  + w_volfrac(i,j ,k+1) * w_arr(i,j ,k+1) + w_volfrac(i,j-1,k+1) * w_arr(i,j-1,k+1))
437  / (w_volfrac(i,j-1,k) + w_volfrac(i,j,k) + w_volfrac(i,j,k+1) + w_volfrac(i,j-1,k+1));
438 
439  // Impose tangential velocity as Dirichlet condition
440  Real v_dot_n = velx * nx + vely * ny + velz * nz;
441  Dirichlet_u = velx - v_dot_n * nx;
442  Dirichlet_v = vely - v_dot_n * ny;
443  Dirichlet_w = velz - v_dot_n * nz;
444  }
445 
446  GpuArray<Real,AMREX_SPACEDIM> slopes_u;
447  GpuArray<Real,AMREX_SPACEDIM> slopes_v;
448  GpuArray<Real,AMREX_SPACEDIM> slopes_w;
449 
450  slopes_u = erf_calc_slopes_eb_Dirichlet_staggered( Vars::yvel, Vars::xvel, dx, dy, dz, i, j, k, bcent_eb, Dirichlet_u, u_arr, u_volcent, u_cellflg);
451  slopes_v = erf_calc_slopes_eb_Dirichlet ( dx, dy, dz, i, j, k, bcent_eb, Dirichlet_v, v_arr, v_volcent, v_cellflg);
452  slopes_w = erf_calc_slopes_eb_Dirichlet_staggered( Vars::yvel, Vars::zvel, dx, dy, dz, i, j, k, bcent_eb, Dirichlet_w, w_arr, w_volcent, w_cellflg);
453 
454  Real dudx = slopes_u[0];
455  Real dudy = slopes_u[1];
456  Real dudz = slopes_u[2];
457  Real dvdx = slopes_v[0];
458  Real dvdy = slopes_v[1];
459  Real dvdz = slopes_v[2];
460  Real dwdx = slopes_w[0];
461  Real dwdy = slopes_w[1];
462  Real dwdz = slopes_w[2];
463 
464  Real tau22_eb = ( dvdy - ( dudx + dvdy + dwdz ) / three );
465  Real tau12_eb = myhalf * (dudy + dvdx);
466  Real tau23_eb = myhalf * (dvdz + dwdy);
467 
468  if (l_no_slip) {
469 
470  dvdn = - mu_eff * (nx * tau12_eb + ny * tau22_eb + nz * tau23_eb);
471 
472  } else if (l_surface_layer) {
473 
474  Real tbx_x, tbx_y, tbx_z, tby_x, tby_y, tby_z;
475  compute_tangent_vectors(nx, ny, nz, tbx_x, tbx_y, tbx_z, tby_x, tby_y, tby_z);
476 
477  Real tau11_eb = ( dudx - ( dudx + dvdy + dwdz ) / three );
478  Real tau33_eb = ( dwdz - ( dudx + dvdy + dwdz ) / three );
479  Real tau13_eb = myhalf * (dudz + dwdx);
480 
481  Real tauzz = mu_eff * ( nx*nx*tau11_eb + ny*ny*tau22_eb + nz*nz*tau33_eb
482  + two * (nx*ny*tau12_eb + ny*nz*tau23_eb + nx*nz*tau13_eb ));
483 
484  dvdn = - tbx_y * v_tau_eb13(i,j,k) - tby_y * v_tau_eb23(i,j,k) - ny * tauzz;
485  }
486  }
487 
488  rho_v_rhs(i,j,k) -= barea * dvdn / (vol * v_volfrac(i,j,k));
489  }
490  }
491  });
492  } // v_type
493 
494  // z-momentum
495  if (w_type == FabType::regular) {
496 
497  ParallelFor(bxz, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
498 
499  Real mfsq = mf_mx(i,j,0) * mf_my(i,j,0);
500 
501  Real diffContrib = ( (tau13(i+1, j , k ) - tau13(i , j , k ) ) * dxinv * mfsq
502  + (tau23(i , j+1, k ) - tau23(i , j , k ) ) * dyinv * mfsq
503  + (tau33(i , j , k ) - tau33(i , j , k-1) ) * dzinv );
504  diffContrib /= w_volfrac(i,j,k);
505 
506  rho_w_rhs(i,j,k) -= diffContrib;
507  });
508 
509  } else if (w_type == FabType::singlevalued) {
510 
511  ParallelFor(bxz, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
512 
513  if (w_volfrac(i,j,k)>zero) {
514 
515  // Inv Jacobian
516  Real mfsq = mf_mx(i,j,0) * mf_my(i,j,0);
517 
518  Real diffContrib = ( (tau13(i+1, j , k ) * w_afrac_x(i+1,j ,k )
519  - tau13(i , j , k ) * w_afrac_x(i ,j ,k ) ) * dxinv * mfsq
520  + (tau23(i , j+1, k ) * w_afrac_y(i ,j+1,k )
521  - tau23(i , j , k ) * w_afrac_y(i ,j ,k ) ) * dyinv * mfsq
522  + (tau33(i , j , k ) * w_afrac_z(i ,j ,k+1)
523  - tau33(i , j , k-1) * w_afrac_z(i ,j ,k ) ) * dzinv );
524  diffContrib /= w_volfrac(i,j,k);
525 
526  rho_w_rhs(i,j,k) -= diffContrib;
527 
528  if (!l_constraint_z && w_cellflg(i,j,k).isSingleValued()) {
529 
530  Real axm = w_afrac_x(i ,j ,k );
531  Real axp = w_afrac_x(i+1,j ,k );
532  Real aym = w_afrac_y(i ,j ,k );
533  Real ayp = w_afrac_y(i ,j+1,k );
534  Real azm = w_afrac_z(i ,j ,k );
535  Real azp = w_afrac_z(i ,j ,k+1);
536 
537  Real adx = (axm-axp) * dy * dz;
538  Real ady = (aym-ayp) * dx * dz;
539  Real adz = (azm-azp) * dx * dy;
540 
541  Real barea = std::sqrt(adx*adx + ady*ady + adz*adz);
542 
543  Real dwdn = zero;
544 
545  if (l_no_slip || l_surface_layer) {
546 
547  const RealVect bcent_eb {w_bcent(i,j,k,0), w_bcent(i,j,k,1), w_bcent(i,j,k,2)};
548 
549  Real Dirichlet_u {zero};
550  Real Dirichlet_v {zero};
551  Real Dirichlet_w {zero};
552 
553  Real nx = w_bnorm(i,j,k,0);
554  Real ny = w_bnorm(i,j,k,1);
555  Real nz = w_bnorm(i,j,k,2);
556 
557  if (l_surface_layer) {
558 
559  // Average u and v onto the z-face
560  Real velx = (u_volfrac(i ,j,k-1) * u_arr(i ,j,k-1) + u_volfrac(i+1,j,k-1) * u_arr(i+1,j,k-1)
561  + u_volfrac(i+1,j,k ) * u_arr(i+1,j,k ) + u_volfrac(i ,j,k ) * u_arr(i ,j,k ))
562  / (u_volfrac(i,j,k-1) + u_volfrac(i+1,j,k-1) + u_volfrac(i+1,j,k) + u_volfrac(i,j,k));
563  Real vely = (v_volfrac(i,j ,k-1) * v_arr(i,j ,k-1) + v_volfrac(i,j+1,k-1) * v_arr(i,j+1,k-1)
564  + v_volfrac(i,j+1,k ) * v_arr(i,j+1,k ) + v_volfrac(i,j ,k ) * v_arr(i,j ,k ))
565  / (v_volfrac(i,j,k-1) + v_volfrac(i,j+1,k-1) + v_volfrac(i,j+1,k) + v_volfrac(i,j,k));
566  Real velz = w_arr(i,j,k);
567 
568  // Impose tangential velocity as Dirichlet condition
569  Real v_dot_n = velx * nx + vely * ny + velz * nz;
570  Dirichlet_u = velx - v_dot_n * nx;
571  Dirichlet_v = vely - v_dot_n * ny;
572  Dirichlet_w = velz - v_dot_n * nz;
573  }
574 
575  GpuArray<Real,AMREX_SPACEDIM> slopes_u;
576  GpuArray<Real,AMREX_SPACEDIM> slopes_v;
577  GpuArray<Real,AMREX_SPACEDIM> slopes_w;
578 
579  slopes_u = erf_calc_slopes_eb_Dirichlet_staggered( Vars::zvel, Vars::xvel, dx, dy, dz, i, j, k, bcent_eb, Dirichlet_u, u_arr, u_volcent, u_cellflg);
580  slopes_v = erf_calc_slopes_eb_Dirichlet_staggered( Vars::zvel, Vars::yvel, dx, dy, dz, i, j, k, bcent_eb, Dirichlet_v, v_arr, v_volcent, v_cellflg);
581  slopes_w = erf_calc_slopes_eb_Dirichlet ( dx, dy, dz, i, j, k, bcent_eb, Dirichlet_w, w_arr, w_volcent, w_cellflg);
582 
583  Real dudx = slopes_u[0];
584  Real dudy = slopes_u[1];
585  Real dudz = slopes_u[2];
586  Real dvdx = slopes_v[0];
587  Real dvdy = slopes_v[1];
588  Real dvdz = slopes_v[2];
589  Real dwdx = slopes_w[0];
590  Real dwdy = slopes_w[1];
591  Real dwdz = slopes_w[2];
592 
593  Real tau33_eb = ( dwdz - ( dudx + dvdy + dwdz ) / three );
594  Real tau13_eb = myhalf * (dudz + dwdx);
595  Real tau23_eb = myhalf * (dvdz + dwdy);
596 
597  if (l_no_slip) {
598 
599  dwdn = - mu_eff * (nx * tau13_eb + ny * tau23_eb + nz * tau33_eb);
600 
601  } else if (l_surface_layer) {
602 
603  Real tbx_x, tbx_y, tbx_z, tby_x, tby_y, tby_z;
604  compute_tangent_vectors(nx, ny, nz, tbx_x, tbx_y, tbx_z, tby_x, tby_y, tby_z);
605 
606  Real tau11_eb = ( dudx - ( dudx + dvdy + dwdz ) / three );
607  Real tau22_eb = ( dvdy - ( dudx + dvdy + dwdz ) / three );
608  Real tau12_eb = myhalf * (dudy + dvdx);
609 
610  Real tauzz = mu_eff * ( nx*nx*tau11_eb + ny*ny*tau22_eb + nz*nz*tau33_eb
611  + two * (nx*ny*tau12_eb + ny*nz*tau23_eb + nx*nz*tau13_eb ));
612 
613  dwdn = - tbx_z * w_tau_eb13(i,j,k) - tby_z * w_tau_eb23(i,j,k) - nz * tauzz;
614 
615  }
616  }
617 
618  rho_w_rhs(i,j,k) -= barea * dwdn / (vol * w_volfrac(i,j,k));
619  }
620  }
621  });
622  } // w_type
623 }
tau12
@ tau12
Definition: ERF_DataStruct.H:40
tau23
@ tau23
Definition: ERF_DataStruct.H:40
tau33
@ tau33
Definition: ERF_DataStruct.H:40
tau22
@ tau22
Definition: ERF_DataStruct.H:40
tau11
@ tau11
Definition: ERF_DataStruct.H:40
tau13
@ tau13
Definition: ERF_DataStruct.H:40
MolecDiffType::ConstantAlpha
@ ConstantAlpha
DiffusionSrcForMom_EB
void DiffusionSrcForMom_EB(const MFIter &mfi, [[maybe_unused]] const Box &domain, const Box &bxx, const Box &bxy, const Box &bxz, const Array4< Real > &rho_u_rhs, const Array4< Real > &rho_v_rhs, const Array4< Real > &rho_w_rhs, const Array4< const Real > &u_arr, const Array4< const Real > &v_arr, const Array4< const Real > &w_arr, const Array4< const Real > &tau11, const Array4< const Real > &tau22, const Array4< const Real > &tau33, const Array4< const Real > &tau12, const Array4< const Real > &tau13, const Array4< const Real > &tau23, const Array4< const Real > &u_tau_eb13, const Array4< const Real > &u_tau_eb23, const Array4< const Real > &v_tau_eb13, const Array4< const Real > &v_tau_eb23, const Array4< const Real > &w_tau_eb13, const Array4< const Real > &w_tau_eb23, const Real *dx_arr, 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 SolverChoice &solverChoice, const eb_ &ebfact, [[maybe_unused]] const BCRec *d_bcrec_ptr)
Definition: ERF_DiffusionSrcForMom_EB.cpp:123
compute_tangent_vectors
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void compute_tangent_vectors(Real nx, Real ny, Real nz, Real &tbx_x, Real &tbx_y, Real &tbx_z, Real &tby_x, Real &tby_y, Real &tby_z)
Definition: ERF_DiffusionSrcForMom_EB.cpp:43
erf_calc_slopes_eb_Dirichlet
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > erf_calc_slopes_eb_Dirichlet(amrex::Real dx, amrex::Real dy, amrex::Real dz, int i, int j, int k, amrex::RealVect const &bcent_eb, amrex::Real const state_eb, amrex::Array4< amrex::Real const > const &state, amrex::Array4< amrex::Real const > const &ccent, amrex::Array4< amrex::EBCellFlag const > const &flag)
Compute least-squares slopes using EB Dirichlet data.
Definition: ERF_EBSlopes.H:28
erf_calc_slopes_eb_Dirichlet_staggered
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > erf_calc_slopes_eb_Dirichlet_staggered(int igrid_query, int igrid_data, amrex::Real dx, amrex::Real dy, amrex::Real dz, int i, int j, int k, amrex::RealVect const &bcent_eb, amrex::Real const state_eb, amrex::Array4< amrex::Real const > const &state, amrex::Array4< amrex::Real const > const &ccent, amrex::Array4< amrex::EBCellFlag const > const &flag)
Compute least-squares slopes from staggered data using EB Dirichlet data.
Definition: ERF_EBSlopes.H:160
EBBoundaryType::NoSlipWall
@ NoSlipWall
EBBoundaryType::SurfaceLayer
@ SurfaceLayer
dxInv
amrex::GpuArray< Real, AMREX_SPACEDIM > dxInv
Definition: ERF_InitCustomPertVels_ParticleTests.H:17
dy
const Real dy
Definition: ERF_InitCustomPert_ABL.H:45
dx
const Real dx
Definition: ERF_InitCustomPert_ABL.H:44
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:32
two
constexpr amrex::Real two
Definition: ERF_NumericalConstants.H:31
myhalf
constexpr amrex::Real myhalf
Definition: ERF_NumericalConstants.H:34
eb_::get_w_const_factory
eb_aux_ const * get_w_const_factory() const noexcept
Return the ERF auxiliary z-face EB factory.
Definition: ERF_EB.H:133
eb_::get_v_const_factory
eb_aux_ const * get_v_const_factory() const noexcept
Return the ERF auxiliary y-face EB factory.
Definition: ERF_EB.H:131
eb_::get_u_const_factory
eb_aux_ const * get_u_const_factory() const noexcept
Return the ERF auxiliary x-face EB factory.
Definition: ERF_EB.H:129
eb_aux_::getMultiEBCellFlagFab
const amrex::FabArray< amrex::EBCellFlagFab > & getMultiEBCellFlagFab() const
Return the reconstructed EB cell flags.
Definition: ERF_EBAux.cpp:1146
Vars::xvel
@ xvel
Definition: ERF_IndexDefines.H:215
Vars::zvel
@ zvel
Definition: ERF_IndexDefines.H:217
Vars::yvel
@ yvel
Definition: ERF_IndexDefines.H:216
WDM6SedCellScratch::dz
@ dz
Definition: ERF_AdvanceWDM6.cpp:272
DiffChoice
Definition: ERF_DiffStruct.H:22
DiffChoice::eb_diff_constraint_z
bool eb_diff_constraint_z
Whether to constrain EB diffusion in the z direction.
Definition: ERF_DiffStruct.H:109
DiffChoice::molec_diff_type
MolecDiffType molec_diff_type
Selected molecular transport model.
Definition: ERF_DiffStruct.H:94
DiffChoice::eb_diff_constraint_y
bool eb_diff_constraint_y
Whether to constrain EB diffusion in the y direction.
Definition: ERF_DiffStruct.H:108
DiffChoice::eb_diff_constraint_x
bool eb_diff_constraint_x
Whether to constrain EB diffusion in the x direction.
Definition: ERF_DiffStruct.H:107
DiffChoice::dynamic_viscosity
amrex::Real dynamic_viscosity
Dynamic viscosity for momentum diffusion [kg/(m-s)].
Definition: ERF_DiffStruct.H:106
DiffChoice::rho0_trans
amrex::Real rho0_trans
Reference density used to compute dynamic diffusion coefficients [kg/m3].
Definition: ERF_DiffStruct.H:101
EBChoice
Definition: ERF_EBStruct.H:36
SolverChoice::ebChoice
EBChoice ebChoice
Embedded-boundary options.
Definition: ERF_DataStruct.H:1975
SolverChoice::diffChoice
DiffChoice diffChoice
Diffusion-related options.
Definition: ERF_DataStruct.H:1971
Here is the call graph for this function: