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