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  constexpr Real tol = Real(1.e-12);
51 
52  // x-tangential vector: t_bx = (e_x - (e_x · n)n) / ||e_x - (e_x · n)n||
53  // e_x = (1,0,0), so e_x · n = nx
54  Real tbx_norm2 = ny*ny + nz*nz;
55  if (tbx_norm2 > tol) {
56  Real tbx_norm_inv = one / std::sqrt(tbx_norm2);
57  tbx_x = ( one - nx * nx) * tbx_norm_inv;
58  tbx_y = ( - nx * ny) * tbx_norm_inv;
59  tbx_z = ( - nx * nz) * tbx_norm_inv;
60  } else {
61  tbx_x = zero;
62  tbx_y = zero;
63  tbx_z = zero;
64  }
65 
66  // y-tangential vector: t_by = (e_y - (e_y · n)n) / ||e_y - (e_y · n)n||
67  // e_y = (0,1,0), so e_y · n = ny
68  Real tby_norm2 = nx*nx + nz*nz;
69  if (tby_norm2 > tol) {
70  Real tby_norm_inv = one / std::sqrt(tby_norm2);
71  tby_x = ( - ny * nx) * tby_norm_inv;
72  tby_y = ( one - ny * ny) * tby_norm_inv;
73  tby_z = ( - ny * nz) * tby_norm_inv;
74  } else {
75  tby_x = zero;
76  tby_y = zero;
77  tby_z = zero;
78  }
79 }
one
constexpr amrex::Real one
Definition: ERF_Constants.H:9
zero
constexpr amrex::Real zero
Definition: ERF_Constants.H:8
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
152 {
153  BL_PROFILE_VAR("DiffusionSrcForMom_EB()",DiffusionSrcForMom_EB);
154 
155  DiffChoice dc = solverChoice.diffChoice;
156  const bool l_use_constAlpha = ( dc.molec_diff_type == MolecDiffType::ConstantAlpha );
157  Real mu_eff = (l_use_constAlpha) ? two * dc.dynamic_viscosity / dc.rho0_trans
158  : two * dc.dynamic_viscosity;
159 
160  auto dxinv = dxInv[0], dyinv = dxInv[1], dzinv = dxInv[2];
161  Real dx = dx_arr[0], dy = dx_arr[1], dz = dx_arr[2];
162  Real vol = dx * dy * dz;
163 
164  EBChoice ebChoice = solverChoice.ebChoice;
165  const bool l_no_slip = (ebChoice.eb_boundary_type == EBBoundaryType::NoSlipWall);
166  const bool l_surface_layer = (ebChoice.eb_boundary_type == EBBoundaryType::SurfaceLayer);
167  const bool l_constraint_x = solverChoice.diffChoice.eb_diff_constraint_x;
168  const bool l_constraint_y = solverChoice.diffChoice.eb_diff_constraint_y;
169  const bool l_constraint_z = solverChoice.diffChoice.eb_diff_constraint_z;
170 
171  // EB u-factory
172  const auto* u_factory = ebfact.get_u_const_factory();
173  Array4<const EBCellFlag> u_cellflg = u_factory->getMultiEBCellFlagFab()[mfi].const_array();
174  Array4<const Real > u_volfrac = u_factory->getVolFrac().const_array(mfi);
175  Array4<const Real > u_volcent{};
176  Array4<const Real > u_afrac_x{};
177  Array4<const Real > u_afrac_y{};
178  Array4<const Real > u_afrac_z{};
179  Array4<const Real > u_bcent{};
180  Array4<const Real > u_bnorm{};
181  FabType u_type = u_factory->getMultiEBCellFlagFab()[mfi].getType();
182  if (u_type == FabType::singlevalued) {
183  u_volcent = u_factory->getCentroid().const_array(mfi);
184  u_afrac_x = u_factory->getAreaFrac()[0]->const_array(mfi);
185  u_afrac_y = u_factory->getAreaFrac()[1]->const_array(mfi);
186  u_afrac_z = u_factory->getAreaFrac()[2]->const_array(mfi);
187  u_bcent = u_factory->getBndryCent().const_array(mfi);
188  u_bnorm = u_factory->getBndryNormal().const_array(mfi);
189  }
190 
191 
192  // EB v-factory
193  const auto* v_factory = ebfact.get_v_const_factory();
194  Array4<const EBCellFlag> v_cellflg = v_factory->getMultiEBCellFlagFab()[mfi].const_array();
195  Array4<const Real > v_volfrac = v_factory->getVolFrac().const_array(mfi);
196  Array4<const Real > v_volcent{};
197  Array4<const Real > v_afrac_x{};
198  Array4<const Real > v_afrac_y{};
199  Array4<const Real > v_afrac_z{};
200  Array4<const Real > v_bcent{};
201  Array4<const Real > v_bnorm{};
202  FabType v_type = v_factory->getMultiEBCellFlagFab()[mfi].getType();
203  if (v_type == FabType::singlevalued) {
204  v_volcent = v_factory->getCentroid().const_array(mfi);
205  v_afrac_x = v_factory->getAreaFrac()[0]->const_array(mfi);
206  v_afrac_y = v_factory->getAreaFrac()[1]->const_array(mfi);
207  v_afrac_z = v_factory->getAreaFrac()[2]->const_array(mfi);
208  v_bcent = v_factory->getBndryCent().const_array(mfi);
209  v_bnorm = v_factory->getBndryNormal().const_array(mfi);
210  }
211 
212  // EB w-factory
213  const auto* w_factory = ebfact.get_w_const_factory();
214  Array4<const EBCellFlag> w_cellflg = w_factory->getMultiEBCellFlagFab()[mfi].const_array();
215  Array4<const Real > w_volfrac = w_factory->getVolFrac().const_array(mfi);
216  Array4<const Real > w_volcent{};
217  Array4<const Real > w_afrac_x{};
218  Array4<const Real > w_afrac_y{};
219  Array4<const Real > w_afrac_z{};
220  Array4<const Real > w_bcent{};
221  Array4<const Real > w_bnorm{};
222  FabType w_type = w_factory->getMultiEBCellFlagFab()[mfi].getType();
223  if (w_type == FabType::singlevalued) {
224  w_volcent = w_factory->getCentroid().const_array(mfi);
225  w_afrac_x = w_factory->getAreaFrac()[0]->const_array(mfi);
226  w_afrac_y = w_factory->getAreaFrac()[1]->const_array(mfi);
227  w_afrac_z = w_factory->getAreaFrac()[2]->const_array(mfi);
228  w_bcent = w_factory->getBndryCent().const_array(mfi);
229  w_bnorm = w_factory->getBndryNormal().const_array(mfi);
230  }
231 
232  // x-momentum
233  if (u_type == FabType::regular) {
234 
235  ParallelFor(bxx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
236 
237  Real mfsq = mf_ux(i,j,0) * mf_uy(i,j,0);
238 
239  Real diffContrib = ( (tau11(i , j , k ) - tau11(i-1, j , k ) ) * dxinv * mfsq
240  + (tau12(i , j+1, k ) - tau12(i , j , k ) ) * dyinv * mfsq
241  + (tau13(i , j , k+1) - tau13(i , j , k ) ) * dzinv );
242  diffContrib /= u_volfrac(i,j,k);
243 
244  rho_u_rhs(i,j,k) -= diffContrib;
245  });
246 
247  } else if (u_type == FabType::singlevalued) {
248 
249  ParallelFor(bxx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
250 
251  if (u_volfrac(i,j,k)>zero) {
252 
253  // Inv Jacobian
254  Real mfsq = mf_ux(i,j,0) * mf_uy(i,j,0);
255 
256  Real diffContrib = ( (tau11(i , j , k ) * u_afrac_x(i+1,j ,k )
257  - tau11(i-1, j , k ) * u_afrac_x(i ,j ,k ) ) * dxinv * mfsq
258  + (tau12(i , j+1, k ) * u_afrac_y(i ,j+1,k )
259  - tau12(i , j , k ) * u_afrac_y(i ,j ,k ) ) * dyinv * mfsq
260  + (tau13(i , j , k+1) * u_afrac_z(i ,j ,k+1)
261  - tau13(i , j , k ) * u_afrac_z(i ,j ,k )) * dzinv );
262  diffContrib /= u_volfrac(i,j,k);
263 
264  rho_u_rhs(i,j,k) -= diffContrib;
265 
266  if (!l_constraint_x && u_cellflg(i,j,k).isSingleValued()) {
267 
268  Real axm = u_afrac_x(i ,j ,k );
269  Real axp = u_afrac_x(i+1,j ,k );
270  Real aym = u_afrac_y(i ,j ,k );
271  Real ayp = u_afrac_y(i ,j+1,k );
272  Real azm = u_afrac_z(i ,j ,k );
273  Real azp = u_afrac_z(i ,j ,k+1);
274 
275  Real adx = (axm-axp) * dy * dz;
276  Real ady = (aym-ayp) * dx * dz;
277  Real adz = (azm-azp) * dx * dy;
278 
279  Real barea = std::sqrt(adx*adx + ady*ady + adz*adz);
280 
281  Real dudn = zero;
282 
283  if (l_no_slip || l_surface_layer) {
284 
285  RealVect bcent_eb {u_bcent(i,j,k,0), u_bcent(i,j,k,1), u_bcent(i,j,k,2)};
286 
287  Real Dirichlet_u {zero};
288  Real Dirichlet_v {zero};
289  Real Dirichlet_w {zero};
290 
291  Real nx = u_bnorm(i,j,k,0);
292  Real ny = u_bnorm(i,j,k,1);
293  Real nz = u_bnorm(i,j,k,2);
294 
295  if (l_surface_layer) {
296 
297  // Average v and w onto the x-face
298  Real velx = u_arr(i,j,k);
299  Real vely = (v_volfrac(i-1,j ,k) * v_arr(i-1,j ,k) + v_volfrac(i,j ,k) * v_arr(i,j ,k)
300  + 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))
301  / (v_volfrac(i-1,j,k) + v_volfrac(i,j,k) + v_volfrac(i-1,j+1,k) + v_volfrac(i,j+1,k));
302 
303  Real velz = (w_volfrac(i-1,j,k ) * w_arr(i-1,j,k ) + w_volfrac(i,j,k ) * w_arr(i,j,k )
304  + 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))
305  / (w_volfrac(i-1,j,k) + w_volfrac(i,j,k) + w_volfrac(i-1,j,k+1) + w_volfrac(i,j,k+1));
306 
307  // Impose tangential velocity as Dirichlet condition
308  Real v_dot_n = velx * nx + vely * ny + velz * nz;
309  Dirichlet_u = velx - v_dot_n * nx;
310  Dirichlet_v = vely - v_dot_n * ny;
311  Dirichlet_w = velz - v_dot_n * nz;
312  }
313 
314  GpuArray<Real,AMREX_SPACEDIM> slopes_u;
315  GpuArray<Real,AMREX_SPACEDIM> slopes_v;
316  GpuArray<Real,AMREX_SPACEDIM> slopes_w;
317 
318  slopes_u = erf_calc_slopes_eb_Dirichlet ( dx, dy, dz, i, j, k, bcent_eb, Dirichlet_u, u_arr, u_volcent, u_cellflg);
319  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);
320  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);
321 
322  Real dudx = slopes_u[0];
323  Real dudy = slopes_u[1];
324  Real dudz = slopes_u[2];
325  Real dvdx = slopes_v[0];
326  Real dvdy = slopes_v[1];
327  Real dvdz = slopes_v[2];
328  Real dwdx = slopes_w[0];
329  Real dwdy = slopes_w[1];
330  Real dwdz = slopes_w[2];
331 
332  Real tau11_eb = ( dudx - ( dudx + dvdy + dwdz ) / three );
333  Real tau12_eb = myhalf * (dudy + dvdx);
334  Real tau13_eb = myhalf * (dudz + dwdx);
335 
336  if (l_no_slip) {
337 
338  dudn = - mu_eff * (nx * tau11_eb + ny * tau12_eb + nz * tau13_eb);
339 
340  } else if (l_surface_layer) {
341 
342  Real tbx_x, tbx_y, tbx_z, tby_x, tby_y, tby_z;
343  compute_tangent_vectors(nx, ny, nz, tbx_x, tbx_y, tbx_z, tby_x, tby_y, tby_z);
344 
345  Real tau22_eb = ( dvdy - ( dudx + dvdy + dwdz ) / three );
346  Real tau33_eb = ( dwdz - ( dudx + dvdy + dwdz ) / three );
347  Real tau23_eb = myhalf * (dvdz + dwdy);
348 
349  Real tauzz = mu_eff * ( nx*nx*tau11_eb + ny*ny*tau22_eb + nz*nz*tau33_eb
350  + two * (nx*ny*tau12_eb + ny*nz*tau23_eb + nx*nz*tau13_eb ));
351 
352  dudn = - tbx_x * u_tau_eb13(i,j,k) - tby_x * u_tau_eb23(i,j,k) - nx * tauzz;
353  }
354  }
355 
356  rho_u_rhs(i,j,k) -= barea * dudn / (vol * u_volfrac(i,j,k));
357  }
358  }
359 
360  });
361  } // u_type
362 
363  // y-momentum
364  if (v_type == FabType::regular) {
365 
366  ParallelFor(bxy, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
367 
368  Real mfsq = mf_vx(i,j,0) * mf_vy(i,j,0);
369 
370  Real diffContrib = ( (tau12(i+1, j , k ) - tau12(i , j , k ) ) * dxinv * mfsq
371  + (tau22(i , j , k ) - tau22(i , j-1, k ) ) * dyinv * mfsq
372  + (tau23(i , j , k+1) - tau23(i , j , k ) ) * dzinv );
373  diffContrib /= v_volfrac(i,j,k);
374 
375  rho_v_rhs(i,j,k) -= diffContrib;
376  });
377  } else if (v_type == FabType::singlevalued) {
378 
379  ParallelFor(bxy, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
380 
381  if (v_volfrac(i,j,k)>zero) {
382 
383  // Inv Jacobian
384  Real mfsq = mf_vx(i,j,0) * mf_vy(i,j,0);
385 
386  Real diffContrib = ( (tau12(i+1, j , k ) * v_afrac_x(i+1,j ,k )
387  - tau12(i , j , k ) * v_afrac_x(i ,j ,k ) ) * dxinv * mfsq
388  + (tau22(i , j , k ) * v_afrac_y(i ,j+1,k )
389  - tau22(i , j-1, k ) * v_afrac_y(i ,j ,k ) ) * dyinv * mfsq
390  + (tau23(i , j , k+1) * v_afrac_z(i ,j ,k+1)
391  - tau23(i , j , k ) * v_afrac_z(i ,j ,k ) ) * dzinv );
392  diffContrib /= v_volfrac(i,j,k);
393 
394  rho_v_rhs(i,j,k) -= diffContrib;
395 
396  if (!l_constraint_y && v_cellflg(i,j,k).isSingleValued()) {
397 
398  Real axm = v_afrac_x(i ,j ,k );
399  Real axp = v_afrac_x(i+1,j ,k );
400  Real aym = v_afrac_y(i ,j ,k );
401  Real ayp = v_afrac_y(i ,j+1,k );
402  Real azm = v_afrac_z(i ,j ,k );
403  Real azp = v_afrac_z(i ,j ,k+1);
404 
405  Real adx = (axm-axp) * dy * dz;
406  Real ady = (aym-ayp) * dx * dz;
407  Real adz = (azm-azp) * dx * dy;
408 
409  Real barea = std::sqrt(adx*adx + ady*ady + adz*adz);
410 
411  Real dvdn = 0.0;
412 
413  if (l_no_slip || l_surface_layer) {
414 
415  RealVect bcent_eb {v_bcent(i,j,k,0), v_bcent(i,j,k,1), v_bcent(i,j,k,2)};
416 
417  Real Dirichlet_u {zero};
418  Real Dirichlet_v {zero};
419  Real Dirichlet_w {zero};
420 
421  Real nx = v_bnorm(i,j,k,0);
422  Real ny = v_bnorm(i,j,k,1);
423  Real nz = v_bnorm(i,j,k,2);
424 
425  if (l_surface_layer) {
426 
427  // Average u and w onto the y-face
428  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)
429  + u_volfrac(i+1,j ,k) * u_arr(i+1,j ,k) + u_volfrac(i ,j ,k) * u_arr(i ,j ,k))
430  / (u_volfrac(i,j-1,k) + u_volfrac(i+1,j-1,k) + u_volfrac(i+1,j,k) + u_volfrac(i,j,k));
431  Real vely = v_arr(i,j,k);
432  Real velz = (w_volfrac(i,j-1,k ) * w_arr(i,j-1,k ) + w_volfrac(i,j,k ) * w_arr(i,j,k )
433  + 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))
434  / (w_volfrac(i,j-1,k) + w_volfrac(i,j,k) + w_volfrac(i,j,k+1) + w_volfrac(i,j-1,k+1));
435 
436  // Impose tangential velocity as Dirichlet condition
437  Real v_dot_n = velx * nx + vely * ny + velz * nz;
438  Dirichlet_u = velx - v_dot_n * nx;
439  Dirichlet_v = vely - v_dot_n * ny;
440  Dirichlet_w = velz - v_dot_n * nz;
441  }
442 
443  GpuArray<Real,AMREX_SPACEDIM> slopes_u;
444  GpuArray<Real,AMREX_SPACEDIM> slopes_v;
445  GpuArray<Real,AMREX_SPACEDIM> slopes_w;
446 
447  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);
448  slopes_v = erf_calc_slopes_eb_Dirichlet ( dx, dy, dz, i, j, k, bcent_eb, Dirichlet_v, v_arr, v_volcent, v_cellflg);
449  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);
450 
451  Real dudx = slopes_u[0];
452  Real dudy = slopes_u[1];
453  Real dudz = slopes_u[2];
454  Real dvdx = slopes_v[0];
455  Real dvdy = slopes_v[1];
456  Real dvdz = slopes_v[2];
457  Real dwdx = slopes_w[0];
458  Real dwdy = slopes_w[1];
459  Real dwdz = slopes_w[2];
460 
461  Real tau22_eb = ( dvdy - ( dudx + dvdy + dwdz ) / three );
462  Real tau12_eb = myhalf * (dudy + dvdx);
463  Real tau23_eb = myhalf * (dvdz + dwdy);
464 
465  if (l_no_slip) {
466 
467  dvdn = - mu_eff * (nx * tau12_eb + ny * tau22_eb + nz * tau23_eb);
468 
469  } else if (l_surface_layer) {
470 
471  Real tbx_x, tbx_y, tbx_z, tby_x, tby_y, tby_z;
472  compute_tangent_vectors(nx, ny, nz, tbx_x, tbx_y, tbx_z, tby_x, tby_y, tby_z);
473 
474  Real tau11_eb = ( dudx - ( dudx + dvdy + dwdz ) / three );
475  Real tau33_eb = ( dwdz - ( dudx + dvdy + dwdz ) / three );
476  Real tau13_eb = myhalf * (dudz + dwdx);
477 
478  Real tauzz = mu_eff * ( nx*nx*tau11_eb + ny*ny*tau22_eb + nz*nz*tau33_eb
479  + two * (nx*ny*tau12_eb + ny*nz*tau23_eb + nx*nz*tau13_eb ));
480 
481  dvdn = - tbx_y * v_tau_eb13(i,j,k) - tby_y * v_tau_eb23(i,j,k) - ny * tauzz;
482  }
483  }
484 
485  rho_v_rhs(i,j,k) -= barea * dvdn / (vol * v_volfrac(i,j,k));
486  }
487  }
488  });
489  } // v_type
490 
491  // z-momentum
492  if (w_type == FabType::regular) {
493 
494  ParallelFor(bxz, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
495 
496  Real mfsq = mf_mx(i,j,0) * mf_my(i,j,0);
497 
498  Real diffContrib = ( (tau13(i+1, j , k ) - tau13(i , j , k ) ) * dxinv * mfsq
499  + (tau23(i , j+1, k ) - tau23(i , j , k ) ) * dyinv * mfsq
500  + (tau33(i , j , k ) - tau33(i , j , k-1) ) * dzinv );
501  diffContrib /= w_volfrac(i,j,k);
502 
503  rho_w_rhs(i,j,k) -= diffContrib;
504  });
505 
506  } else if (w_type == FabType::singlevalued) {
507 
508  ParallelFor(bxz, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
509 
510  if (w_volfrac(i,j,k)>zero) {
511 
512  // Inv Jacobian
513  Real mfsq = mf_mx(i,j,0) * mf_my(i,j,0);
514 
515  Real diffContrib = ( (tau13(i+1, j , k ) * w_afrac_x(i+1,j ,k )
516  - tau13(i , j , k ) * w_afrac_x(i ,j ,k ) ) * dxinv * mfsq
517  + (tau23(i , j+1, k ) * w_afrac_y(i ,j+1,k )
518  - tau23(i , j , k ) * w_afrac_y(i ,j ,k ) ) * dyinv * mfsq
519  + (tau33(i , j , k ) * w_afrac_z(i ,j ,k+1)
520  - tau33(i , j , k-1) * w_afrac_z(i ,j ,k ) ) * dzinv );
521  diffContrib /= w_volfrac(i,j,k);
522 
523  rho_w_rhs(i,j,k) -= diffContrib;
524 
525  if (!l_constraint_z && w_cellflg(i,j,k).isSingleValued()) {
526 
527  Real axm = w_afrac_x(i ,j ,k );
528  Real axp = w_afrac_x(i+1,j ,k );
529  Real aym = w_afrac_y(i ,j ,k );
530  Real ayp = w_afrac_y(i ,j+1,k );
531  Real azm = w_afrac_z(i ,j ,k );
532  Real azp = w_afrac_z(i ,j ,k+1);
533 
534  Real adx = (axm-axp) * dy * dz;
535  Real ady = (aym-ayp) * dx * dz;
536  Real adz = (azm-azp) * dx * dy;
537 
538  Real barea = std::sqrt(adx*adx + ady*ady + adz*adz);
539 
540  Real dwdn = zero;
541 
542  if (l_no_slip || l_surface_layer) {
543 
544  const RealVect bcent_eb {w_bcent(i,j,k,0), w_bcent(i,j,k,1), w_bcent(i,j,k,2)};
545 
546  Real Dirichlet_u {zero};
547  Real Dirichlet_v {zero};
548  Real Dirichlet_w {zero};
549 
550  Real nx = w_bnorm(i,j,k,0);
551  Real ny = w_bnorm(i,j,k,1);
552  Real nz = w_bnorm(i,j,k,2);
553 
554  if (l_surface_layer) {
555 
556  // Average u and v onto the z-face
557  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)
558  + u_volfrac(i+1,j,k ) * u_arr(i+1,j,k ) + u_volfrac(i ,j,k ) * u_arr(i ,j,k ))
559  / (u_volfrac(i,j,k-1) + u_volfrac(i+1,j,k-1) + u_volfrac(i+1,j,k) + u_volfrac(i,j,k));
560  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)
561  + v_volfrac(i,j+1,k ) * v_arr(i,j+1,k ) + v_volfrac(i,j ,k ) * v_arr(i,j ,k ))
562  / (v_volfrac(i,j,k-1) + v_volfrac(i,j+1,k-1) + v_volfrac(i,j+1,k) + v_volfrac(i,j,k));
563  Real velz = w_arr(i,j,k);
564 
565  // Impose tangential velocity as Dirichlet condition
566  Real v_dot_n = velx * nx + vely * ny + velz * nz;
567  Dirichlet_u = velx - v_dot_n * nx;
568  Dirichlet_v = vely - v_dot_n * ny;
569  Dirichlet_w = velz - v_dot_n * nz;
570  }
571 
572  GpuArray<Real,AMREX_SPACEDIM> slopes_u;
573  GpuArray<Real,AMREX_SPACEDIM> slopes_v;
574  GpuArray<Real,AMREX_SPACEDIM> slopes_w;
575 
576  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);
577  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);
578  slopes_w = erf_calc_slopes_eb_Dirichlet ( dx, dy, dz, i, j, k, bcent_eb, Dirichlet_w, w_arr, w_volcent, w_cellflg);
579 
580  Real dudx = slopes_u[0];
581  Real dudy = slopes_u[1];
582  Real dudz = slopes_u[2];
583  Real dvdx = slopes_v[0];
584  Real dvdy = slopes_v[1];
585  Real dvdz = slopes_v[2];
586  Real dwdx = slopes_w[0];
587  Real dwdy = slopes_w[1];
588  Real dwdz = slopes_w[2];
589 
590  Real tau33_eb = ( dwdz - ( dudx + dvdy + dwdz ) / three );
591  Real tau13_eb = myhalf * (dudz + dwdx);
592  Real tau23_eb = myhalf * (dvdz + dwdy);
593 
594  if (l_no_slip) {
595 
596  dwdn = - mu_eff * (nx * tau13_eb + ny * tau23_eb + nz * tau33_eb);
597 
598  } else if (l_surface_layer) {
599 
600  Real tbx_x, tbx_y, tbx_z, tby_x, tby_y, tby_z;
601  compute_tangent_vectors(nx, ny, nz, tbx_x, tbx_y, tbx_z, tby_x, tby_y, tby_z);
602 
603  Real tau11_eb = ( dudx - ( dudx + dvdy + dwdz ) / three );
604  Real tau22_eb = ( dvdy - ( dudx + dvdy + dwdz ) / three );
605  Real tau12_eb = myhalf * (dudy + dvdx);
606 
607  Real tauzz = mu_eff * ( nx*nx*tau11_eb + ny*ny*tau22_eb + nz*nz*tau33_eb
608  + two * (nx*ny*tau12_eb + ny*nz*tau23_eb + nx*nz*tau13_eb ));
609 
610  dwdn = - tbx_z * w_tau_eb13(i,j,k) - tby_z * w_tau_eb23(i,j,k) - nz * tauzz;
611 
612  }
613  }
614 
615  rho_w_rhs(i,j,k) -= barea * dwdn / (vol * w_volfrac(i,j,k));
616  }
617  }
618  });
619  } // w_type
620 }
three
constexpr amrex::Real three
Definition: ERF_Constants.H:11
two
constexpr amrex::Real two
Definition: ERF_Constants.H:10
myhalf
constexpr amrex::Real myhalf
Definition: ERF_Constants.H:13
tau12
@ tau12
Definition: ERF_DataStruct.H:39
tau23
@ tau23
Definition: ERF_DataStruct.H:39
tau33
@ tau33
Definition: ERF_DataStruct.H:39
tau22
@ tau22
Definition: ERF_DataStruct.H:39
tau11
@ tau11
Definition: ERF_DataStruct.H:39
tau13
@ tau13
Definition: ERF_DataStruct.H:39
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:120
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);})
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:123
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:121
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:119
eb_aux_::getMultiEBCellFlagFab
const amrex::FabArray< amrex::EBCellFlagFab > & getMultiEBCellFlagFab() const
Return the reconstructed EB cell flags.
Definition: ERF_EBAux.cpp:1145
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:270
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:1865
SolverChoice::diffChoice
DiffChoice diffChoice
Diffusion-related options.
Definition: ERF_DataStruct.H:1861
Here is the call graph for this function: