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 orthonormal tangent vectors at an EB boundary given the normal vector. Uses Gram-Schmidt orthogonalization against standard basis vectors.

Parameters
[in]nx,ny,nzComponents of the normal vector
[out]tbx_x,tbx_y,tbx_zComponents of first tangent vector (from e_x)
[out]tby_x,tby_y,tby_zComponents of second tangent vector (from e_y)
24 {
25  // x-tangential vector: t_bx = (e_x - (e_x · n)n) / ||e_x - (e_x · n)n||
26  // e_x = (1,0,0), so e_x · n = nx
27  tbx_x = one - nx * nx;
28  tbx_y = - nx * ny;
29  tbx_z = - nx * nz;
30  Real tbx_norm = std::sqrt(tbx_x*tbx_x + tbx_y*tbx_y + tbx_z*tbx_z);
31  tbx_x /= tbx_norm;
32  tbx_y /= tbx_norm;
33  tbx_z /= tbx_norm;
34 
35  // y-tangential vector: t_by = (e_y - (e_y · n)n) / ||e_y - (e_y · n)n||
36  // e_y = (0,1,0), so e_y · n = ny
37  tby_x = - ny * nx;
38  tby_y = one - ny * ny;
39  tby_z = - ny * nz;
40  Real tby_norm = std::sqrt(tby_x*tby_x + tby_y*tby_y + tby_z*tby_z);
41  tby_x /= tby_norm;
42  tby_y /= tby_norm;
43  tby_z /= tby_norm;
44 }
one
constexpr amrex::Real one
Definition: ERF_Constants.H:7
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]tau1111 stress
[in]tau2222 stress
[in]tau3333 stress
[in]tau1212 stress
[in]tau1313 stress
[in]tau2323 stress
[in]dxInvinverse cell size array
[in]mf_mmap factor at cell center
[in]ebfactEB factories for cell- and face-centered variables
100 {
101  BL_PROFILE_VAR("DiffusionSrcForMom_EB()",DiffusionSrcForMom_EB);
102 
103  DiffChoice dc = solverChoice.diffChoice;
104  const bool l_use_constAlpha = ( dc.molec_diff_type == MolecDiffType::ConstantAlpha );
105  Real mu_eff = (l_use_constAlpha) ? two * dc.dynamic_viscosity / dc.rho0_trans
106  : two * dc.dynamic_viscosity;
107 
108  auto dxinv = dxInv[0], dyinv = dxInv[1], dzinv = dxInv[2];
109  Real dx = dx_arr[0], dy = dx_arr[1], dz = dx_arr[2];
110  Real vol = dx * dy * dz;
111 
112  EBChoice ebChoice = solverChoice.ebChoice;
113  const bool l_no_slip = (ebChoice.eb_boundary_type == EBBoundaryType::NoSlipWall);
114  const bool l_surface_layer = (ebChoice.eb_boundary_type == EBBoundaryType::SurfaceLayer);
115  const bool l_constraint_x = solverChoice.diffChoice.eb_diff_constraint_x;
116  const bool l_constraint_y = solverChoice.diffChoice.eb_diff_constraint_y;
117  const bool l_constraint_z = solverChoice.diffChoice.eb_diff_constraint_z;
118 
119  // int n = 0;
120 
121  // EB u-factory
122  const auto* u_factory = ebfact.get_u_const_factory();
123  Array4<const EBCellFlag> u_cellflg = u_factory->getMultiEBCellFlagFab()[mfi].const_array();
124  Array4<const Real > u_volfrac = u_factory->getVolFrac().const_array(mfi);
125  Array4<const Real > u_volcent = u_factory->getCentroid().const_array(mfi);
126  Array4<const Real > u_afrac_x = u_factory->getAreaFrac()[0]->const_array(mfi);
127  Array4<const Real > u_afrac_y = u_factory->getAreaFrac()[1]->const_array(mfi);
128  Array4<const Real > u_afrac_z = u_factory->getAreaFrac()[2]->const_array(mfi);
129  Array4<const Real > u_bcent = u_factory->getBndryCent().const_array(mfi);
130  Array4<const Real > u_bnorm = u_factory->getBndryNorm().const_array(mfi);
131 
132  // EB v-factory
133  const auto* v_factory = ebfact.get_v_const_factory();
134  Array4<const EBCellFlag> v_cellflg = v_factory->getMultiEBCellFlagFab()[mfi].const_array();
135  Array4<const Real > v_volfrac = v_factory->getVolFrac().const_array(mfi);
136  Array4<const Real > v_volcent = v_factory->getCentroid().const_array(mfi);
137  Array4<const Real > v_afrac_x = v_factory->getAreaFrac()[0]->const_array(mfi);
138  Array4<const Real > v_afrac_y = v_factory->getAreaFrac()[1]->const_array(mfi);
139  Array4<const Real > v_afrac_z = v_factory->getAreaFrac()[2]->const_array(mfi);
140  Array4<const Real > v_bcent = v_factory->getBndryCent().const_array(mfi);
141  Array4<const Real > v_bnorm = v_factory->getBndryNorm().const_array(mfi);
142 
143  // EB w-factory
144  const auto* w_factory = ebfact.get_w_const_factory();
145  Array4<const EBCellFlag> w_cellflg = w_factory->getMultiEBCellFlagFab()[mfi].const_array();
146  Array4<const Real > w_volfrac = w_factory->getVolFrac().const_array(mfi);
147  Array4<const Real > w_volcent = w_factory->getCentroid().const_array(mfi);
148  Array4<const Real > w_afrac_x = w_factory->getAreaFrac()[0]->const_array(mfi);
149  Array4<const Real > w_afrac_y = w_factory->getAreaFrac()[1]->const_array(mfi);
150  Array4<const Real > w_afrac_z = w_factory->getAreaFrac()[2]->const_array(mfi);
151  Array4<const Real > w_bcent = w_factory->getBndryCent().const_array(mfi);
152  Array4<const Real > w_bnorm = w_factory->getBndryNorm().const_array(mfi);
153 
154  // x-momentum
155  ParallelFor(bxx,
156  [=] AMREX_GPU_DEVICE (int i, int j, int k)
157  {
158  if (u_volfrac(i,j,k)>zero) {
159 
160  // Inv Jacobian
161  Real mfsq = mf_ux(i,j,0) * mf_uy(i,j,0);
162 
163  Real diffContrib = ( (tau11(i , j , k ) * u_afrac_x(i+1,j ,k )
164  - tau11(i-1, j , k ) * u_afrac_x(i ,j ,k ) ) * dxinv * mfsq
165  + (tau12(i , j+1, k ) * u_afrac_y(i ,j+1,k )
166  - tau12(i , j , k ) * u_afrac_y(i ,j ,k ) ) * dyinv * mfsq
167  + (tau13(i , j , k+1) * u_afrac_z(i ,j ,k+1)
168  - tau13(i , j , k ) * u_afrac_z(i ,j ,k )) * dzinv );
169  diffContrib /= u_volfrac(i,j,k);
170 
171  rho_u_rhs(i,j,k) -= diffContrib;
172 
173  if (!l_constraint_x && u_cellflg(i,j,k).isSingleValued()) {
174 
175  Real axm = u_afrac_x(i ,j ,k );
176  Real axp = u_afrac_x(i+1,j ,k );
177  Real aym = u_afrac_y(i ,j ,k );
178  Real ayp = u_afrac_y(i ,j+1,k );
179  Real azm = u_afrac_z(i ,j ,k );
180  Real azp = u_afrac_z(i ,j ,k+1);
181 
182  Real adx = (axm-axp) * dy * dz;
183  Real ady = (aym-ayp) * dx * dz;
184  Real adz = (azm-azp) * dx * dy;
185 
186  Real barea = std::sqrt(adx*adx + ady*ady + adz*adz);
187 
188  Real dudn = zero;
189 
190  if (l_no_slip || l_surface_layer) {
191 
192  RealVect bcent_eb {u_bcent(i,j,k,0), u_bcent(i,j,k,1), u_bcent(i,j,k,2)};
193 
194  Real Dirichlet_u {zero};
195  Real Dirichlet_v {zero};
196  Real Dirichlet_w {zero};
197 
198  Real nx = u_bnorm(i,j,k,0);
199  Real ny = u_bnorm(i,j,k,1);
200  Real nz = u_bnorm(i,j,k,2);
201 
202  if (l_surface_layer) {
203 
204  // Average v and w onto the x-face
205  Real velx = u_arr(i,j,k);
206  Real vely = (v_volfrac(i-1,j ,k) * v_arr(i-1,j ,k) + v_volfrac(i,j ,k) * v_arr(i,j ,k)
207  + 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))
208  / (v_volfrac(i-1,j,k) + v_volfrac(i,j,k) + v_volfrac(i-1,j+1,k) + v_volfrac(i,j+1,k));
209 
210  Real velz = (w_volfrac(i-1,j,k ) * w_arr(i-1,j,k ) + w_volfrac(i,j,k ) * w_arr(i,j,k )
211  + 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))
212  / (w_volfrac(i-1,j,k) + w_volfrac(i,j,k) + w_volfrac(i-1,j,k+1) + w_volfrac(i,j,k+1));
213 
214  // Impose tangential velocity as Dirichlet condition
215  Real v_dot_n = velx * nx + vely * ny + velz * nz;
216  Dirichlet_u = velx - v_dot_n * nx;
217  Dirichlet_v = vely - v_dot_n * ny;
218  Dirichlet_w = velz - v_dot_n * nz;
219  }
220 
221  GpuArray<Real,AMREX_SPACEDIM> slopes_u;
222  GpuArray<Real,AMREX_SPACEDIM> slopes_v;
223  GpuArray<Real,AMREX_SPACEDIM> slopes_w;
224 
225  slopes_u = erf_calc_slopes_eb_Dirichlet ( dx, dy, dz, i, j, k, bcent_eb, Dirichlet_u, u_arr, u_volcent, u_cellflg);
226  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);
227  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);
228 
229  Real dudx = slopes_u[0];
230  Real dudy = slopes_u[1];
231  Real dudz = slopes_u[2];
232  Real dvdx = slopes_v[0];
233  Real dvdy = slopes_v[1];
234  Real dvdz = slopes_v[2];
235  Real dwdx = slopes_w[0];
236  Real dwdy = slopes_w[1];
237  Real dwdz = slopes_w[2];
238 
239  Real tau11_eb = ( dudx - ( dudx + dvdy + dwdz ) / three );
240  Real tau12_eb = myhalf * (dudy + dvdx);
241  Real tau13_eb = myhalf * (dudz + dwdx);
242 
243  if (l_no_slip) {
244 
245  dudn = - mu_eff * (nx * tau11_eb + ny * tau12_eb + nz * tau13_eb);
246 
247  } else if (l_surface_layer) {
248 
249  Real tbx_x, tbx_y, tbx_z, tby_x, tby_y, tby_z;
250  compute_tangent_vectors(nx, ny, nz, tbx_x, tbx_y, tbx_z, tby_x, tby_y, tby_z);
251 
252  Real tau22_eb = ( dvdy - ( dudx + dvdy + dwdz ) / three );
253  Real tau33_eb = ( dwdz - ( dudx + dvdy + dwdz ) / three );
254  Real tau23_eb = myhalf * (dvdz + dwdy);
255 
256  Real tauzz = mu_eff * ( nx*nx*tau11_eb + ny*ny*tau22_eb + nz*nz*tau33_eb
257  + two * (nx*ny*tau12_eb + ny*nz*tau23_eb + nx*nz*tau13_eb ));
258 
259  dudn = - tbx_x * u_tau_eb13(i,j,k) - tby_x * u_tau_eb23(i,j,k) - nx * tauzz;
260  }
261  }
262 
263  rho_u_rhs(i,j,k) -= barea * dudn / (vol * u_volfrac(i,j,k));
264  }
265  }
266 
267  });
268 
269  // y-momentum
270  ParallelFor(bxy,
271  [=] AMREX_GPU_DEVICE (int i, int j, int k)
272  {
273  if (v_volfrac(i,j,k)>zero) {
274 
275  // Inv Jacobian
276  Real mfsq = mf_vx(i,j,0) * mf_vy(i,j,0);
277 
278  Real diffContrib = ( (tau12(i+1, j , k ) * v_afrac_x(i+1,j ,k )
279  - tau12(i , j , k ) * v_afrac_x(i ,j ,k ) ) * dxinv * mfsq
280  + (tau22(i , j , k ) * v_afrac_y(i ,j+1,k )
281  - tau22(i , j-1, k ) * v_afrac_y(i ,j ,k ) ) * dyinv * mfsq
282  + (tau23(i , j , k+1) * v_afrac_z(i ,j ,k+1)
283  - tau23(i , j , k ) * v_afrac_z(i ,j ,k ) ) * dzinv );
284  diffContrib /= v_volfrac(i,j,k);
285 
286  rho_v_rhs(i,j,k) -= diffContrib;
287 
288  if (!l_constraint_y && v_cellflg(i,j,k).isSingleValued()) {
289 
290  Real axm = v_afrac_x(i ,j ,k );
291  Real axp = v_afrac_x(i+1,j ,k );
292  Real aym = v_afrac_y(i ,j ,k );
293  Real ayp = v_afrac_y(i ,j+1,k );
294  Real azm = v_afrac_z(i ,j ,k );
295  Real azp = v_afrac_z(i ,j ,k+1);
296 
297  Real adx = (axm-axp) * dy * dz;
298  Real ady = (aym-ayp) * dx * dz;
299  Real adz = (azm-azp) * dx * dy;
300 
301  Real barea = std::sqrt(adx*adx + ady*ady + adz*adz);
302 
303  Real dvdn = 0.0;
304 
305  if (l_no_slip || l_surface_layer) {
306 
307  RealVect bcent_eb {v_bcent(i,j,k,0), v_bcent(i,j,k,1), v_bcent(i,j,k,2)};
308 
309  Real Dirichlet_u {zero};
310  Real Dirichlet_v {zero};
311  Real Dirichlet_w {zero};
312 
313  Real nx = v_bnorm(i,j,k,0);
314  Real ny = v_bnorm(i,j,k,1);
315  Real nz = v_bnorm(i,j,k,2);
316 
317  if (l_surface_layer) {
318 
319  // Average u and w onto the y-face
320  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)
321  + u_volfrac(i+1,j ,k) * u_arr(i+1,j ,k) + u_volfrac(i ,j ,k) * u_arr(i ,j ,k))
322  / (u_volfrac(i,j-1,k) + u_volfrac(i+1,j-1,k) + u_volfrac(i+1,j,k) + u_volfrac(i,j,k));
323  Real vely = v_arr(i,j,k);
324  Real velz = (w_volfrac(i,j-1,k ) * w_arr(i,j-1,k ) + w_volfrac(i,j,k ) * w_arr(i,j,k )
325  + 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))
326  / (w_volfrac(i,j-1,k) + w_volfrac(i,j,k) + w_volfrac(i,j,k+1) + w_volfrac(i,j-1,k+1));
327 
328  // Impose tangential velocity as Dirichlet condition
329  Real v_dot_n = velx * nx + vely * ny + velz * nz;
330  Dirichlet_u = velx - v_dot_n * nx;
331  Dirichlet_v = vely - v_dot_n * ny;
332  Dirichlet_w = velz - v_dot_n * nz;
333  }
334 
335  GpuArray<Real,AMREX_SPACEDIM> slopes_u;
336  GpuArray<Real,AMREX_SPACEDIM> slopes_v;
337  GpuArray<Real,AMREX_SPACEDIM> slopes_w;
338 
339  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);
340  slopes_v = erf_calc_slopes_eb_Dirichlet ( dx, dy, dz, i, j, k, bcent_eb, Dirichlet_v, v_arr, v_volcent, v_cellflg);
341  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);
342 
343  Real dudx = slopes_u[0];
344  Real dudy = slopes_u[1];
345  Real dudz = slopes_u[2];
346  Real dvdx = slopes_v[0];
347  Real dvdy = slopes_v[1];
348  Real dvdz = slopes_v[2];
349  Real dwdx = slopes_w[0];
350  Real dwdy = slopes_w[1];
351  Real dwdz = slopes_w[2];
352 
353  Real tau22_eb = ( dvdy - ( dudx + dvdy + dwdz ) / three );
354  Real tau12_eb = myhalf * (dudy + dvdx);
355  Real tau23_eb = myhalf * (dvdz + dwdy);
356 
357  if (l_no_slip) {
358 
359  dvdn = - mu_eff * (nx * tau12_eb + ny * tau22_eb + nz * tau23_eb);
360 
361  } else if (l_surface_layer) {
362 
363  Real tbx_x, tbx_y, tbx_z, tby_x, tby_y, tby_z;
364  compute_tangent_vectors(nx, ny, nz, tbx_x, tbx_y, tbx_z, tby_x, tby_y, tby_z);
365 
366  Real tau11_eb = ( dudx - ( dudx + dvdy + dwdz ) / three );
367  Real tau33_eb = ( dwdz - ( dudx + dvdy + dwdz ) / three );
368  Real tau13_eb = myhalf * (dudz + dwdx);
369 
370  Real tauzz = mu_eff * ( nx*nx*tau11_eb + ny*ny*tau22_eb + nz*nz*tau33_eb
371  + two * (nx*ny*tau12_eb + ny*nz*tau23_eb + nx*nz*tau13_eb ));
372 
373  dvdn = - tbx_y * v_tau_eb13(i,j,k) - tby_y * v_tau_eb23(i,j,k) - ny * tauzz;
374  }
375  }
376 
377  rho_v_rhs(i,j,k) -= barea * dvdn / (vol * v_volfrac(i,j,k));
378  }
379  }
380  });
381 
382  // z-momentum
383  ParallelFor(bxz,
384  [=] AMREX_GPU_DEVICE (int i, int j, int k)
385  {
386  if (w_volfrac(i,j,k)>zero) {
387 
388  // Inv Jacobian
389  Real mfsq = mf_mx(i,j,0) * mf_my(i,j,0);
390 
391  Real diffContrib = ( (tau13(i+1, j , k ) * w_afrac_x(i+1,j ,k )
392  - tau13(i , j , k ) * w_afrac_x(i ,j ,k ) ) * dxinv * mfsq
393  + (tau23(i , j+1, k ) * w_afrac_y(i ,j+1,k )
394  - tau23(i , j , k ) * w_afrac_y(i ,j ,k ) ) * dyinv * mfsq
395  + (tau33(i , j , k ) * w_afrac_z(i ,j ,k+1)
396  - tau33(i , j , k-1) * w_afrac_z(i ,j ,k ) ) * dzinv );
397  diffContrib /= w_volfrac(i,j,k);
398  rho_w_rhs(i,j,k) -= diffContrib;
399 
400  if (!l_constraint_z && w_cellflg(i,j,k).isSingleValued()) {
401 
402  Real axm = w_afrac_x(i ,j ,k );
403  Real axp = w_afrac_x(i+1,j ,k );
404  Real aym = w_afrac_y(i ,j ,k );
405  Real ayp = w_afrac_y(i ,j+1,k );
406  Real azm = w_afrac_z(i ,j ,k );
407  Real azp = w_afrac_z(i ,j ,k+1);
408 
409  Real adx = (axm-axp) * dy * dz;
410  Real ady = (aym-ayp) * dx * dz;
411  Real adz = (azm-azp) * dx * dy;
412 
413  Real barea = std::sqrt(adx*adx + ady*ady + adz*adz);
414 
415  Real dwdn = zero;
416 
417  if (l_no_slip || l_surface_layer) {
418 
419  const RealVect bcent_eb {w_bcent(i,j,k,0), w_bcent(i,j,k,1), w_bcent(i,j,k,2)};
420 
421  Real Dirichlet_u {zero};
422  Real Dirichlet_v {zero};
423  Real Dirichlet_w {zero};
424 
425  Real nx = w_bnorm(i,j,k,0);
426  Real ny = w_bnorm(i,j,k,1);
427  Real nz = w_bnorm(i,j,k,2);
428 
429  if (l_surface_layer) {
430 
431  // Average u and v onto the z-face
432  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)
433  + u_volfrac(i+1,j,k ) * u_arr(i+1,j,k ) + u_volfrac(i ,j,k ) * u_arr(i ,j,k ))
434  / (u_volfrac(i,j,k-1) + u_volfrac(i+1,j,k-1) + u_volfrac(i+1,j,k) + u_volfrac(i,j,k));
435  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)
436  + v_volfrac(i,j+1,k ) * v_arr(i,j+1,k ) + v_volfrac(i,j ,k ) * v_arr(i,j ,k ))
437  / (v_volfrac(i,j,k-1) + v_volfrac(i,j+1,k-1) + v_volfrac(i,j+1,k) + v_volfrac(i,j,k));
438  Real velz = w_arr(i,j,k);
439 
440  // Impose tangential velocity as Dirichlet condition
441  Real v_dot_n = velx * nx + vely * ny + velz * nz;
442  Dirichlet_u = velx - v_dot_n * nx;
443  Dirichlet_v = vely - v_dot_n * ny;
444  Dirichlet_w = velz - v_dot_n * nz;
445  }
446 
447  GpuArray<Real,AMREX_SPACEDIM> slopes_u;
448  GpuArray<Real,AMREX_SPACEDIM> slopes_v;
449  GpuArray<Real,AMREX_SPACEDIM> slopes_w;
450 
451  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);
452  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);
453  slopes_w = erf_calc_slopes_eb_Dirichlet ( dx, dy, dz, i, j, k, bcent_eb, Dirichlet_w, w_arr, w_volcent, w_cellflg);
454 
455  Real dudx = slopes_u[0];
456  Real dudy = slopes_u[1];
457  Real dudz = slopes_u[2];
458  Real dvdx = slopes_v[0];
459  Real dvdy = slopes_v[1];
460  Real dvdz = slopes_v[2];
461  Real dwdx = slopes_w[0];
462  Real dwdy = slopes_w[1];
463  Real dwdz = slopes_w[2];
464 
465  Real tau33_eb = ( dwdz - ( dudx + dvdy + dwdz ) / three );
466  Real tau13_eb = myhalf * (dudz + dwdx);
467  Real tau23_eb = myhalf * (dvdz + dwdy);
468 
469  if (l_no_slip) {
470 
471  dwdn = - mu_eff * (nx * tau13_eb + ny * tau23_eb + nz * tau33_eb);
472 
473  } else if (l_surface_layer) {
474 
475  Real tbx_x, tbx_y, tbx_z, tby_x, tby_y, tby_z;
476  compute_tangent_vectors(nx, ny, nz, tbx_x, tbx_y, tbx_z, tby_x, tby_y, tby_z);
477 
478  Real tau11_eb = ( dudx - ( dudx + dvdy + dwdz ) / three );
479  Real tau22_eb = ( dvdy - ( dudx + dvdy + dwdz ) / three );
480  Real tau12_eb = myhalf * (dudy + dvdx);
481 
482  Real tauzz = mu_eff * ( nx*nx*tau11_eb + ny*ny*tau22_eb + nz*nz*tau33_eb
483  + two * (nx*ny*tau12_eb + ny*nz*tau23_eb + nx*nz*tau13_eb ));
484 
485  dwdn = - tbx_z * w_tau_eb13(i,j,k) - tby_z * w_tau_eb23(i,j,k) - nz * tauzz;
486 
487  }
488  }
489 
490  rho_w_rhs(i,j,k) -= barea * dwdn / (vol * w_volfrac(i,j,k));
491  }
492  }
493  });
494 
495 }
three
constexpr amrex::Real three
Definition: ERF_Constants.H:9
two
constexpr amrex::Real two
Definition: ERF_Constants.H:8
zero
constexpr amrex::Real zero
Definition: ERF_Constants.H:6
myhalf
constexpr amrex::Real myhalf
Definition: ERF_Constants.H:11
tau12
@ tau12
Definition: ERF_DataStruct.H:32
tau23
@ tau23
Definition: ERF_DataStruct.H:32
tau33
@ tau33
Definition: ERF_DataStruct.H:32
tau22
@ tau22
Definition: ERF_DataStruct.H:32
tau11
@ tau11
Definition: ERF_DataStruct.H:32
tau13
@ tau13
Definition: ERF_DataStruct.H:32
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:68
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:21
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)
Definition: ERF_EBSlopes.H:20
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)
Definition: ERF_EBSlopes.H:146
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:24
dx
const Real dx
Definition: ERF_InitCustomPert_ABL.H:23
ParallelFor
ParallelFor(grown_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
Definition: ERF_EB.H:52
eb_::get_v_const_factory
eb_aux_ const * get_v_const_factory() const noexcept
Definition: ERF_EB.H:51
eb_::get_u_const_factory
eb_aux_ const * get_u_const_factory() const noexcept
Definition: ERF_EB.H:50
eb_aux_::getMultiEBCellFlagFab
const amrex::FabArray< amrex::EBCellFlagFab > & getMultiEBCellFlagFab() const
Definition: ERF_EBAux.cpp:1144
Vars::xvel
@ xvel
Definition: ERF_IndexDefines.H:175
Vars::zvel
@ zvel
Definition: ERF_IndexDefines.H:177
Vars::yvel
@ yvel
Definition: ERF_IndexDefines.H:176
WSM6SedCellScratch::dz
@ dz
Definition: ERF_AdvanceWSM6.cpp:104
DiffChoice
Definition: ERF_DiffStruct.H:19
DiffChoice::rho0_trans
amrex::Real rho0_trans
Definition: ERF_DiffStruct.H:91
DiffChoice::eb_diff_constraint_z
bool eb_diff_constraint_z
Definition: ERF_DiffStruct.H:99
DiffChoice::eb_diff_constraint_x
bool eb_diff_constraint_x
Definition: ERF_DiffStruct.H:97
DiffChoice::eb_diff_constraint_y
bool eb_diff_constraint_y
Definition: ERF_DiffStruct.H:98
DiffChoice::molec_diff_type
MolecDiffType molec_diff_type
Definition: ERF_DiffStruct.H:84
DiffChoice::dynamic_viscosity
amrex::Real dynamic_viscosity
Definition: ERF_DiffStruct.H:96
EBChoice
Definition: ERF_EBStruct.H:27
SolverChoice::diffChoice
DiffChoice diffChoice
Definition: ERF_DataStruct.H:1143
SolverChoice::ebChoice
EBChoice ebChoice
Definition: ERF_DataStruct.H:1147
Here is the call graph for this function: