ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
Functions
ERF_AdvectionSrcForMom_EB.cpp File Reference
#include "AMReX_BCRec.H"
#include <ERF_Advection.H>
#include <ERF_AdvectionSrcForMom_N.H>
#include <ERF_AdvectionSrcForMom_T.H>
#include <ERF_EBAdvectionSrcForMom.H>
Include dependency graph for ERF_AdvectionSrcForMom_EB.cpp:

Functions

void AdvectionSrcForMom_EB (const MFIter &mfi, const Box &bxx, const Box &bxy, const Box &bxz, const Vector< Box > &bxx_grown, const Vector< Box > &bxy_grown, const Vector< Box > &bxz_grown, const Array4< Real > &rho_u_rhs, const Array4< Real > &rho_v_rhs, const Array4< Real > &rho_w_rhs, const Array4< const Real > &u, const Array4< const Real > &v, const Array4< const Real > &w, const Array4< const Real > &rho_u, const Array4< const Real > &rho_v, const Array4< const Real > &omega, const GpuArray< Real, AMREX_SPACEDIM > &cellSizeInv, 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 AdvType horiz_adv_type, const AdvType vert_adv_type, const Real horiz_upw_frac, const Real vert_upw_frac, const eb_ &ebfact, GpuArray< Array4< Real >, AMREX_SPACEDIM > &flx_u_arr, GpuArray< Array4< Real >, AMREX_SPACEDIM > &flx_v_arr, GpuArray< Array4< Real >, AMREX_SPACEDIM > &flx_w_arr, const Vector< iMultiFab > &physbnd_mask, const bool already_on_centroids, const int lo_z_face, const int hi_z_face, const Box &)
 

Function Documentation

◆ AdvectionSrcForMom_EB()

void AdvectionSrcForMom_EB ( const MFIter &  mfi,
const Box &  bxx,
const Box &  bxy,
const Box &  bxz,
const Vector< Box > &  bxx_grown,
const Vector< Box > &  bxy_grown,
const Vector< Box > &  bxz_grown,
const Array4< Real > &  rho_u_rhs,
const Array4< Real > &  rho_v_rhs,
const Array4< Real > &  rho_w_rhs,
const Array4< const Real > &  u,
const Array4< const Real > &  v,
const Array4< const Real > &  w,
const Array4< const Real > &  rho_u,
const Array4< const Real > &  rho_v,
const Array4< const Real > &  omega,
const GpuArray< Real, AMREX_SPACEDIM > &  cellSizeInv,
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 AdvType  horiz_adv_type,
const AdvType  vert_adv_type,
const Real  horiz_upw_frac,
const Real  vert_upw_frac,
const eb_ &  ebfact,
GpuArray< Array4< Real >, AMREX_SPACEDIM > &  flx_u_arr,
GpuArray< Array4< Real >, AMREX_SPACEDIM > &  flx_v_arr,
GpuArray< Array4< Real >, AMREX_SPACEDIM > &  flx_w_arr,
const Vector< iMultiFab > &  physbnd_mask,
const bool  already_on_centroids,
const int  lo_z_face,
const int  hi_z_face,
const Box &   
)

Function for computing the advective tendency for the momentum equations when using EB

Parameters
[in]mfiMultiFab Iterator
[in]bxxbox over which the x-momentum is updated
[in]bxybox over which the y-momentum is updated
[in]bxzbox over which the z-momentum is updated
[in]bxx_growngrown boxes of bxx to loop over the nodal grids of bxx
[in]bxy_growngrown boxes of bxy to loop over the nodal grids of bxy
[in]bxz_growngrown boxes of bxz to loop over the nodal grids of bxz
[out]rho_u_rhstendency for the x-momentum equation
[out]rho_v_rhstendency for the y-momentum equation
[out]rho_w_rhstendency for the z-momentum equation
[in]ux-component of the velocity
[in]vy-component of the velocity
[in]wz-component of the velocity
[in]rho_ux-component of the momentum
[in]rho_vy-component of the momentum
[in]omegacomponent of the momentum normal to the z-coordinate surface
[in]cellSizeInvinverse of the grid spacing
[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]horiz_adv_typesets the spatial order to be used for lateral derivatives
[in]vert_adv_typesets the spatial order to be used for vertical derivatives
[in]horiz_upw_frachorizontal upwind blending fraction
[in]vert_upw_fracvertical upwind blending fraction
[in]ebfactEB factories for cell- and face-centered variables
[in,out]flx_u_arrunweighted flux densities for x-momentum (zero on closed faces)
[in,out]flx_v_arrunweighted flux densities for y-momentum (zero on closed faces)
[in,out]flx_w_arrunweighted flux densities for z-momentum (zero on closed faces)
[in]physbnd_maskVector of masks for flux interpolation (=1 otherwise, =0 if physbnd)
[in]already_on_centroidsflag whether flux interpolation is unnecessary
[in]lo_z_faceminimum z-face k-index at this level
[in]hi_z_facemaximum z-face k-index at this level
84 {
85  BL_PROFILE_VAR("AdvectionSrcForMom_EB", AdvectionSrcForMom_EB);
86 
87  AMREX_ALWAYS_ASSERT(bxz.smallEnd(2) > 0);
88 
89  auto dxInv = cellSizeInv[0], dyInv = cellSizeInv[1], dzInv = cellSizeInv[2];
90 
91  // compute mapfactor inverses
92  Box box2d_u(bxx); box2d_u.setRange(2,0); box2d_u.grow({3,3,0});
93  Box box2d_v(bxy); box2d_v.setRange(2,0); box2d_v.grow({3,3,0});
94 
95  FArrayBox mf_uy_invFAB(box2d_u,1,The_Async_Arena());
96  FArrayBox mf_vx_invFAB(box2d_v,1,The_Async_Arena());
97  const Array4<Real>& mf_uy_inv = mf_uy_invFAB.array();
98  const Array4<Real>& mf_vx_inv = mf_vx_invFAB.array();
99 
100  ParallelFor(box2d_u, box2d_v,
101  [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept
102  {
103  mf_uy_inv(i,j,0) = one / mf_uy(i,j,0);
104  },
105  [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept
106  {
107  mf_vx_inv(i,j,0) = one / mf_vx(i,j,0);
108  });
109 
110  // EB u-factory
111  auto const* u_factory = ebfact.get_u_const_factory();
112  Array4<const EBCellFlag> u_cflag = u_factory->getMultiEBCellFlagFab()[mfi].const_array();
113  Array4<const Real> u_vfrac = u_factory->getVolFrac().const_array(mfi);
114  Array4<const Real> u_afrac_x{};
115  Array4<const Real> u_afrac_y{};
116  Array4<const Real> u_afrac_z{};
117  Array4<const Real> u_fcx{};
118  Array4<const Real> u_fcy{};
119  Array4<const Real> u_fcz{};
120  FabType u_type = u_factory->getMultiEBCellFlagFab()[mfi].getType(bxx);
121 
122  if (u_type == FabType::singlevalued) {
123  u_afrac_x = u_factory->getAreaFrac()[0]->const_array(mfi);
124  u_afrac_y = u_factory->getAreaFrac()[1]->const_array(mfi);
125  u_afrac_z = u_factory->getAreaFrac()[2]->const_array(mfi);
126  u_fcx = u_factory->getFaceCent()[0]->const_array(mfi);
127  u_fcy = u_factory->getFaceCent()[1]->const_array(mfi);
128  u_fcz = u_factory->getFaceCent()[2]->const_array(mfi);
129  }
130 
131  // EB v-factory
132  auto const* v_factory = ebfact.get_v_const_factory();
133  Array4<const EBCellFlag> v_cflag = v_factory->getMultiEBCellFlagFab()[mfi].const_array();
134  Array4<const Real> v_vfrac = v_factory->getVolFrac().const_array(mfi);
135  Array4<const Real> v_afrac_x{};
136  Array4<const Real> v_afrac_y{};
137  Array4<const Real> v_afrac_z{};
138  Array4<const Real> v_fcx{};
139  Array4<const Real> v_fcy{};
140  Array4<const Real> v_fcz{};
141  FabType v_type = v_factory->getMultiEBCellFlagFab()[mfi].getType();
142  if (v_type == FabType::singlevalued) {
143  v_afrac_x = v_factory->getAreaFrac()[0]->const_array(mfi);
144  v_afrac_y = v_factory->getAreaFrac()[1]->const_array(mfi);
145  v_afrac_z = v_factory->getAreaFrac()[2]->const_array(mfi);
146  v_fcx = v_factory->getFaceCent()[0]->const_array(mfi);
147  v_fcy = v_factory->getFaceCent()[1]->const_array(mfi);
148  v_fcz = v_factory->getFaceCent()[2]->const_array(mfi);
149  }
150 
151  // EB w-factory
152  auto const* w_factory = ebfact.get_w_const_factory();
153  Array4<const EBCellFlag> w_cflag = w_factory->getMultiEBCellFlagFab()[mfi].const_array();
154  Array4<const Real> w_vfrac = w_factory->getVolFrac().const_array(mfi);
155  Array4<const Real> w_afrac_x{};
156  Array4<const Real> w_afrac_y{};
157  Array4<const Real> w_afrac_z{};
158  Array4<const Real> w_fcx{};
159  Array4<const Real> w_fcy{};
160  Array4<const Real> w_fcz{};
161  FabType w_type = w_factory->getMultiEBCellFlagFab()[mfi].getType();
162  if (w_type == FabType::singlevalued) {
163  w_afrac_x = w_factory->getAreaFrac()[0]->const_array(mfi);
164  w_afrac_y = w_factory->getAreaFrac()[1]->const_array(mfi);
165  w_afrac_z = w_factory->getAreaFrac()[2]->const_array(mfi);
166  w_fcx = w_factory->getFaceCent()[0]->const_array(mfi);
167  w_fcy = w_factory->getFaceCent()[1]->const_array(mfi);
168  w_fcz = w_factory->getFaceCent()[2]->const_array(mfi);
169  }
170 
171  // Store unweighted flux densities on open faces. Interpolate to face
172  // centroids before applying area fractions in the divergence below.
173  // Inline with 2nd order for efficiency
174  if (horiz_adv_type == AdvType::Centered_2nd && vert_adv_type == AdvType::Centered_2nd)
175  {
176  // Fluxes for x-momentum
177  ParallelFor(bxx_grown[0], bxx_grown[1], bxx_grown[2],
178  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
179  {
180  if (u_type != FabType::covered) {
181  Real flux_base = fourth * (rho_u(i,j,k) * mf_uy_inv(i,j,0) + rho_u(i-1,j,k) * mf_uy_inv(i-1,j,0))
182  * (u(i-1,j,k) + u(i,j,k));
183  if (u_type == FabType::regular) {
184  flx_u_arr[0](i,j,k) = flux_base;
185  } else { // singlevalued
186  flx_u_arr[0](i,j,k) = (u_afrac_x(i,j,k) > zero) ? flux_base : zero;
187  }
188  } else {
189  flx_u_arr[0](i,j,k) = zero;
190  }
191  },
192  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
193  {
194  if (u_type != FabType::covered) {
195  Real flux_base = fourth * (rho_v(i,j,k) * mf_vx_inv(i,j,0) + rho_v(i-1,j,k) * mf_vx_inv(i-1,j,0))
196  * (u(i,j-1,k) + u(i,j,k));
197  if (u_type == FabType::regular) {
198  flx_u_arr[1](i,j,k) = flux_base;
199  } else { // singlevalued
200  flx_u_arr[1](i,j,k) = (u_afrac_y(i,j,k) > zero) ? flux_base : zero;
201  }
202  } else {
203  flx_u_arr[1](i,j,k) = zero;
204  }
205  },
206  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
207  {
208  if (u_type != FabType::covered) {
209  Real flux_base = fourth * (omega(i,j,k) + omega(i-1,j,k)) * (u(i,j,k-1) + u(i,j,k));
210  if (u_type == FabType::regular) {
211  flx_u_arr[2](i,j,k) = flux_base;
212  } else { // singlevalued
213  flx_u_arr[2](i,j,k) = (u_afrac_z(i,j,k) > zero) ? flux_base : zero;
214  }
215  } else {
216  flx_u_arr[2](i,j,k) = zero;
217  }
218  });
219  // Fluxes for y-momentum
220  ParallelFor(bxy_grown[0], bxy_grown[1], bxy_grown[2],
221  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
222  {
223  if (v_type != FabType::covered) {
224  Real flux_base = fourth * (rho_u(i,j,k) * mf_uy_inv(i,j,0) + rho_u(i,j-1,k) * mf_uy_inv(i,j-1,0))
225  * (v(i-1,j,k) + v(i,j,k));
226  if (v_type == FabType::regular) {
227  flx_v_arr[0](i,j,k) = flux_base;
228  } else { // singlevalued
229  flx_v_arr[0](i,j,k) = (v_afrac_x(i,j,k) > zero) ? flux_base : zero;
230  }
231  } else {
232  flx_v_arr[0](i,j,k) = zero;
233  }
234  },
235  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
236  {
237  if (v_type != FabType::covered) {
238  Real flux_base = fourth * (rho_v(i,j,k) * mf_vx_inv(i,j,0) + rho_v(i,j-1,k) * mf_vx_inv(i,j-1,0))
239  * (v(i,j-1,k) + v(i,j,k));
240  if (v_type == FabType::regular) {
241  flx_v_arr[1](i,j,k) = flux_base;
242  } else { // singlevalued
243  flx_v_arr[1](i,j,k) = (v_afrac_y(i,j,k) > zero) ? flux_base : zero;
244  }
245  } else {
246  flx_v_arr[1](i,j,k) = zero;
247  }
248  },
249  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
250  {
251  if (v_type != FabType::covered) {
252  Real flux_base = fourth * (omega(i,j,k) + omega(i,j-1,k)) * (v(i,j,k-1) + v(i,j,k));
253  if (v_type == FabType::regular) {
254  flx_v_arr[2](i,j,k) = flux_base;
255  } else { // singlevalued
256  flx_v_arr[2](i,j,k) = (v_afrac_z(i,j,k) > zero) ? flux_base : zero;
257  }
258  } else {
259  flx_v_arr[2](i,j,k) = zero;
260  }
261  });
262  // Fluxes for z-momentum
263  ParallelFor(bxz_grown[0], bxz_grown[1], bxz_grown[2],
264  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
265  {
266  if (w_type != FabType::covered) {
267  Real flux_base = fourth * (rho_u(i,j,k) + rho_u(i,j, k-1)) * mf_uy_inv(i,j,0)
268  * (w(i-1,j,k) + w(i,j,k));
269  if (w_type == FabType::regular) {
270  flx_w_arr[0](i,j,k) = flux_base;
271  } else { // singlevalued
272  flx_w_arr[0](i,j,k) = (w_afrac_x(i,j,k) > zero) ? flux_base : zero;
273  }
274  } else {
275  flx_w_arr[0](i,j,k) = zero;
276  }
277  },
278  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
279  {
280  if (w_type != FabType::covered) {
281  Real flux_base = fourth * (rho_v(i,j,k) + rho_v(i,j,k-1)) * mf_vx_inv(i,j,0)
282  * (w(i,j-1,k) + w(i,j,k));
283  if (w_type == FabType::regular) {
284  flx_w_arr[1](i,j,k) = flux_base;
285  } else { // singlevalued
286  flx_w_arr[1](i,j,k) = (w_afrac_y(i,j,k) > zero) ? flux_base : zero;
287  }
288  } else {
289  flx_w_arr[1](i,j,k) = zero;
290  }
291  },
292  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
293  {
294  if (w_type != FabType::covered) {
295  Real flux_base = (k==hi_z_face+1) ? omega(i,j,k) * w(i,j,k) :
296  fourth * (omega(i,j,k) + omega(i,j,k-1)) * (w(i,j,k) + w(i,j,k-1));
297  if (w_type == FabType::regular) {
298  flx_w_arr[2](i,j,k) = flux_base;
299  } else { // singlevalued
300  flx_w_arr[2](i,j,k) = (w_afrac_z(i,j,k) > zero) ? flux_base : zero;
301  }
302  } else {
303  flx_w_arr[2](i,j,k) = zero;
304  }
305  });
306 
307  // Template higher order methods
308  } else {
309 
310  if (horiz_adv_type == AdvType::Centered_2nd) {
311  EBAdvectionSrcForMomVert<CENTERED2>(bxx_grown, bxy_grown, bxz_grown,
312  rho_u, rho_v, omega, u, v, w,
313  u_type, v_type, w_type,
314  u_cflag, u_afrac_x, u_afrac_y, u_afrac_z,
315  v_cflag, v_afrac_x, v_afrac_y, v_afrac_z,
316  w_cflag, w_afrac_x, w_afrac_y, w_afrac_z,
317  mf_uy_inv, mf_vx_inv,
318  horiz_upw_frac, vert_upw_frac, vert_adv_type,
319  flx_u_arr, flx_v_arr, flx_w_arr,
320  lo_z_face, hi_z_face);
321  } else if (horiz_adv_type == AdvType::Upwind_3rd) {
322  EBAdvectionSrcForMomVert<UPWIND3>( bxx_grown, bxy_grown, bxz_grown,
323  rho_u, rho_v, omega, u, v, w,
324  u_type, v_type, w_type,
325  u_cflag, u_afrac_x, u_afrac_y, u_afrac_z,
326  v_cflag, v_afrac_x, v_afrac_y, v_afrac_z,
327  w_cflag, w_afrac_x, w_afrac_y, w_afrac_z,
328  mf_uy_inv, mf_vx_inv,
329  horiz_upw_frac, vert_upw_frac, vert_adv_type,
330  flx_u_arr, flx_v_arr, flx_w_arr,
331  lo_z_face, hi_z_face);
332  } else if (horiz_adv_type == AdvType::Centered_4th) {
333  EBAdvectionSrcForMomVert<CENTERED4>(bxx_grown, bxy_grown, bxz_grown,
334  rho_u, rho_v, omega, u, v, w,
335  u_type, v_type, w_type,
336  u_cflag, u_afrac_x, u_afrac_y, u_afrac_z,
337  v_cflag, v_afrac_x, v_afrac_y, v_afrac_z,
338  w_cflag, w_afrac_x, w_afrac_y, w_afrac_z,
339  mf_uy_inv, mf_vx_inv,
340  horiz_upw_frac, vert_upw_frac, vert_adv_type,
341  flx_u_arr, flx_v_arr, flx_w_arr,
342  lo_z_face, hi_z_face);
343  } else if (horiz_adv_type == AdvType::Upwind_5th) {
344  EBAdvectionSrcForMomVert<UPWIND5>( bxx_grown, bxy_grown, bxz_grown,
345  rho_u, rho_v, omega, u, v, w,
346  u_type, v_type, w_type,
347  u_cflag, u_afrac_x, u_afrac_y, u_afrac_z,
348  v_cflag, v_afrac_x, v_afrac_y, v_afrac_z,
349  w_cflag, w_afrac_x, w_afrac_y, w_afrac_z,
350  mf_uy_inv, mf_vx_inv,
351  horiz_upw_frac, vert_upw_frac, vert_adv_type,
352  flx_u_arr, flx_v_arr, flx_w_arr,
353  lo_z_face, hi_z_face);
354  } else if (horiz_adv_type == AdvType::Centered_6th) {
355  EBAdvectionSrcForMomVert<CENTERED6>(bxx_grown, bxy_grown, bxz_grown,
356  rho_u, rho_v, omega, u, v, w,
357  u_type, v_type, w_type,
358  u_cflag, u_afrac_x, u_afrac_y, u_afrac_z,
359  v_cflag, v_afrac_x, v_afrac_y, v_afrac_z,
360  w_cflag, w_afrac_x, w_afrac_y, w_afrac_z,
361  mf_uy_inv, mf_vx_inv,
362  horiz_upw_frac, vert_upw_frac, vert_adv_type,
363  flx_u_arr, flx_v_arr, flx_w_arr,
364  lo_z_face, hi_z_face);
365  } else {
366  AMREX_ASSERT_WITH_MESSAGE(false, "Unknown advection scheme!");
367  }
368  } // horiz_adv_type
369 
370  // Update momentum RHS using the fluxes
371  if (already_on_centroids) {
372 
373  ParallelFor(bxx,
374  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
375  {
376  if (u_vfrac(i,j,k)>zero) {
377  Real mfsq = mf_ux(i,j,0) * mf_uy(i,j,0);
378 
379  // Regular fabs have unit face areas and no area-fraction arrays.
380  Real axm = (u_type == FabType::singlevalued) ? u_afrac_x(i,j,k) : one;
381  Real axp = (u_type == FabType::singlevalued) ? u_afrac_x(i+1,j,k) : one;
382  Real aym = (u_type == FabType::singlevalued) ? u_afrac_y(i,j,k) : one;
383  Real ayp = (u_type == FabType::singlevalued) ? u_afrac_y(i,j+1,k) : one;
384  Real azm = (u_type == FabType::singlevalued) ? u_afrac_z(i,j,k) : one;
385  Real azp = (u_type == FabType::singlevalued) ? u_afrac_z(i,j,k+1) : one;
386 
387  Real advectionSrc = ( (axp * flx_u_arr[0](i+1, j , k ) - axm * flx_u_arr[0](i, j, k)) * dxInv * mfsq
388  + (ayp * flx_u_arr[1](i , j+1, k ) - aym * flx_u_arr[1](i, j, k)) * dyInv * mfsq
389  + (azp * flx_u_arr[2](i , j , k+1) - azm * flx_u_arr[2](i, j, k)) * dzInv ) / u_vfrac(i,j,k);
390  rho_u_rhs(i, j, k) = -advectionSrc;
391  } else {
392  rho_u_rhs(i, j, k) = zero;
393  }
394  });
395 
396  ParallelFor(bxy,
397  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
398  {
399  if (v_vfrac(i,j,k)>zero) {
400  Real mfsq = mf_vx(i,j,0) * mf_vy(i,j,0);
401 
402  // Regular fabs have unit face areas and no area-fraction arrays.
403  Real axm = (v_type == FabType::singlevalued) ? v_afrac_x(i,j,k) : one;
404  Real axp = (v_type == FabType::singlevalued) ? v_afrac_x(i+1,j,k) : one;
405  Real aym = (v_type == FabType::singlevalued) ? v_afrac_y(i,j,k) : one;
406  Real ayp = (v_type == FabType::singlevalued) ? v_afrac_y(i,j+1,k) : one;
407  Real azm = (v_type == FabType::singlevalued) ? v_afrac_z(i,j,k) : one;
408  Real azp = (v_type == FabType::singlevalued) ? v_afrac_z(i,j,k+1) : one;
409 
410  Real advectionSrc = ( (axp * flx_v_arr[0](i+1, j , k ) - axm * flx_v_arr[0](i, j, k)) * dxInv * mfsq
411  + (ayp * flx_v_arr[1](i , j+1, k ) - aym * flx_v_arr[1](i, j, k)) * dyInv * mfsq
412  + (azp * flx_v_arr[2](i , j , k+1) - azm * flx_v_arr[2](i, j, k)) * dzInv ) / v_vfrac(i,j,k);
413  rho_v_rhs(i, j, k) = -advectionSrc;
414  } else {
415  rho_v_rhs(i, j, k) = zero;
416  }
417  });
418 
419  ParallelFor(bxz,
420  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
421  {
422  if (w_vfrac(i,j,k)>zero) {
423  Real mfsq = mf_mx(i,j,0) * mf_my(i,j,0);
424 
425  // Regular fabs have unit face areas and no area-fraction arrays.
426  Real axm = (w_type == FabType::singlevalued) ? w_afrac_x(i,j,k) : one;
427  Real axp = (w_type == FabType::singlevalued) ? w_afrac_x(i+1,j,k) : one;
428  Real aym = (w_type == FabType::singlevalued) ? w_afrac_y(i,j,k) : one;
429  Real ayp = (w_type == FabType::singlevalued) ? w_afrac_y(i,j+1,k) : one;
430  Real azm = (w_type == FabType::singlevalued) ? w_afrac_z(i,j,k) : one;
431  Real azp = (w_type == FabType::singlevalued) ? w_afrac_z(i,j,k+1) : one;
432 
433  Real advectionSrc = ( (axp * flx_w_arr[0](i+1, j , k ) - axm * flx_w_arr[0](i, j, k)) * dxInv * mfsq
434  + (ayp * flx_w_arr[1](i , j+1, k ) - aym * flx_w_arr[1](i, j, k)) * dyInv * mfsq
435  + (azp * flx_w_arr[2](i , j , k+1) - azm * flx_w_arr[2](i, j, k)) * dzInv ) / w_vfrac(i,j,k);
436  rho_w_rhs(i, j, k) = -advectionSrc;
437  } else {
438  rho_w_rhs(i, j, k) = zero;
439  }
440  });
441 
442  } else {
443  // !already_on_centroids
444 
445  Array4<const int> u_mask = physbnd_mask[IntVars::xmom].const_array(mfi);
446  Array4<const int> v_mask = physbnd_mask[IntVars::ymom].const_array(mfi);
447  Array4<const int> w_mask = physbnd_mask[IntVars::zmom].const_array(mfi);
448 
449  if (u_type == FabType::covered) {
450 
451  ParallelFor(bxx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
452  rho_u_rhs(i, j, k) = zero;
453  });
454 
455  } else if (u_type == FabType::regular) {
456 
457  ParallelFor(bxx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
458  Real mfsq = mf_ux(i,j,0) * mf_uy(i,j,0);
459  rho_u_rhs(i, j, k) = - ( (flx_u_arr[0](i+1, j , k ) - flx_u_arr[0](i, j, k)) * dxInv * mfsq
460  + (flx_u_arr[1](i , j+1, k ) - flx_u_arr[1](i, j, k)) * dyInv * mfsq
461  + (flx_u_arr[2](i , j , k+1) - flx_u_arr[2](i, j, k)) * dzInv );
462  });
463 
464  } else if (u_type == FabType::singlevalued) {
465 
466  ParallelFor(bxx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
467 
468  if (u_vfrac(i,j,k)>zero) {
469  Real mfsq = mf_ux(i,j,0) * mf_uy(i,j,0);
470 
471  if (u_cflag(i,j,k).isCovered()) {
472 
473  rho_u_rhs(i, j, k) = zero;
474 
475  } else if (u_cflag(i,j,k).isRegular()) {
476 
477  rho_u_rhs(i, j, k) = - ( (u_afrac_x(i+1, j , k ) * flx_u_arr[0](i+1, j , k ) - u_afrac_x(i, j, k) * flx_u_arr[0](i, j, k)) * dxInv * mfsq
478  + (u_afrac_y(i , j+1, k ) * flx_u_arr[1](i , j+1, k ) - u_afrac_y(i, j, k) * flx_u_arr[1](i, j, k)) * dyInv * mfsq
479  + (u_afrac_z(i , j , k+1) * flx_u_arr[2](i , j , k+1) - u_afrac_z(i, j, k) * flx_u_arr[2](i, j, k)) * dzInv ) / u_vfrac(i,j,k);
480  } else {
481 
482  // Bilinear interpolation
483  Real fxm = flx_u_arr[0](i,j,k);
484  if (u_afrac_x(i,j,k) != zero && u_afrac_x(i,j,k) != one) {
485  int jj = j + static_cast<int>(std::copysign(one, u_fcx(i,j,k,0)));
486  int kk = k + static_cast<int>(std::copysign(one, u_fcx(i,j,k,1)));
487  Real fracy = (u_mask(i-1,jj,k) || u_mask(i,jj,k)) ? std::abs(u_fcx(i,j,k,0)) : zero;
488  Real fracz = (u_mask(i-1,j,kk) || u_mask(i,j,kk)) ? std::abs(u_fcx(i,j,k,1)) : zero;
489  fxm = (one-fracy)*(one-fracz)*fxm
490  + fracy *(one-fracz)*flx_u_arr[0](i,jj,k )
491  + fracz *(one-fracy)*flx_u_arr[0](i,j ,kk)
492  + fracy * fracz *flx_u_arr[0](i,jj,kk);
493  }
494 
495  Real fxp = flx_u_arr[0](i+1,j,k);
496  if (u_afrac_x(i+1,j,k) != zero && u_afrac_x(i+1,j,k) != one) {
497  int jj = j + static_cast<int>(std::copysign(one,u_fcx(i+1,j,k,0)));
498  int kk = k + static_cast<int>(std::copysign(one,u_fcx(i+1,j,k,1)));
499  Real fracy = (u_mask(i,jj,k) || u_mask(i+1,jj,k)) ? std::abs(u_fcx(i+1,j,k,0)) : zero;
500  Real fracz = (u_mask(i,j,kk) || u_mask(i+1,j,kk)) ? std::abs(u_fcx(i+1,j,k,1)) : zero;
501  fxp = (one-fracy)*(one-fracz)*fxp
502  + fracy *(one-fracz)*flx_u_arr[0](i+1,jj,k )
503  + fracz *(one-fracy)*flx_u_arr[0](i+1,j ,kk)
504  + fracy * fracz *flx_u_arr[0](i+1,jj,kk);
505  }
506 
507  Real fym = flx_u_arr[1](i,j,k);
508  if (u_afrac_y(i,j,k) != zero && u_afrac_y(i,j,k) != one) {
509  int ii = i + static_cast<int>(std::copysign(one,u_fcy(i,j,k,0)));
510  int kk = k + static_cast<int>(std::copysign(one,u_fcy(i,j,k,1)));
511  Real fracx = (u_mask(ii,j-1,k) || u_mask(ii,j,k)) ? std::abs(u_fcy(i,j,k,0)) : zero;
512  Real fracz = (u_mask(i,j-1,kk) || u_mask(i,j,kk)) ? std::abs(u_fcy(i,j,k,1)) : zero;
513  fym = (one-fracx)*(one-fracz)*fym
514  + fracx *(one-fracz)*flx_u_arr[1](ii,j,k )
515  + fracz *(one-fracx)*flx_u_arr[1](i ,j,kk)
516  + fracx * fracz *flx_u_arr[1](ii,j,kk);
517  }
518 
519  Real fyp = flx_u_arr[1](i,j+1,k);
520  if (u_afrac_y(i,j+1,k) != zero && u_afrac_y(i,j+1,k) != one) {
521  int ii = i + static_cast<int>(std::copysign(one,u_fcy(i,j+1,k,0)));
522  int kk = k + static_cast<int>(std::copysign(one,u_fcy(i,j+1,k,1)));
523  Real fracx = (u_mask(ii,j,k) || u_mask(ii,j+1,k)) ? std::abs(u_fcy(i,j+1,k,0)) : zero;
524  Real fracz = (u_mask(i,j,kk) || u_mask(i,j+1,kk)) ? std::abs(u_fcy(i,j+1,k,1)) : zero;
525  fyp = (one-fracx)*(one-fracz)*fyp
526  + fracx *(one-fracz)*flx_u_arr[1](ii,j+1,k )
527  + fracz *(one-fracx)*flx_u_arr[1](i ,j+1,kk)
528  + fracx * fracz *flx_u_arr[1](ii,j+1,kk);
529  }
530 
531  Real fzm = flx_u_arr[2](i,j,k);
532  if (u_afrac_z(i,j,k) != zero && u_afrac_z(i,j,k) != one) {
533  int ii = i + static_cast<int>(std::copysign(one,u_fcz(i,j,k,0)));
534  int jj = j + static_cast<int>(std::copysign(one,u_fcz(i,j,k,1)));
535  Real fracx = (u_mask(ii,j,k-1) || u_mask(ii,j,k)) ? std::abs(u_fcz(i,j,k,0)) : zero;
536  Real fracy = (u_mask(i,jj,k-1) || u_mask(i,jj,k)) ? std::abs(u_fcz(i,j,k,1)) : zero;
537  fzm = (one-fracx)*(one-fracy)*fzm
538  + fracx *(one-fracy)*flx_u_arr[2](ii,j ,k)
539  + fracy *(one-fracx)*flx_u_arr[2](i ,jj,k)
540  + fracx * fracy *flx_u_arr[2](ii,jj,k);
541  }
542 
543  Real fzp = flx_u_arr[2](i,j,k+1);
544  if (u_afrac_z(i,j,k+1) != zero && u_afrac_z(i,j,k+1) != one) {
545  int ii = i + static_cast<int>(std::copysign(one,u_fcz(i,j,k+1,0)));
546  int jj = j + static_cast<int>(std::copysign(one,u_fcz(i,j,k+1,1)));
547  Real fracx = (u_mask(ii,j,k) || u_mask(ii,j,k+1)) ? std::abs(u_fcz(i,j,k+1,0)) : zero;
548  Real fracy = (u_mask(i,jj,k) || u_mask(i,jj,k+1)) ? std::abs(u_fcz(i,j,k+1,1)) : zero;
549  fzp = (one-fracx)*(one-fracy)*fzp
550  + fracx *(one-fracy)*flx_u_arr[2](ii,j ,k+1)
551  + fracy *(one-fracx)*flx_u_arr[2](i ,jj,k+1)
552  + fracx * fracy *flx_u_arr[2](ii,jj,k+1);
553  }
554 
555  rho_u_rhs(i, j, k) = - ( (u_afrac_x(i+1, j , k ) * fxp - u_afrac_x(i, j, k) * fxm) * dxInv * mfsq
556  + (u_afrac_y(i , j+1, k ) * fyp - u_afrac_y(i, j, k) * fym) * dyInv * mfsq
557  + (u_afrac_z(i , j , k+1) * fzp - u_afrac_z(i, j, k) * fzm) * dzInv ) / u_vfrac(i,j,k);
558  }
559 
560  } else {
561  rho_u_rhs(i, j, k) = zero;
562  }
563  });
564  } // u_type
565 
566  if (v_type == FabType::covered) {
567 
568  ParallelFor(bxy, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
569  rho_v_rhs(i, j, k) = zero;
570  });
571 
572  } else if (v_type == FabType::regular) {
573 
574  ParallelFor(bxy, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
575  Real mfsq = mf_vx(i,j,0) * mf_vy(i,j,0);
576  rho_v_rhs(i, j, k) = - ( (flx_v_arr[0](i+1, j , k ) - flx_v_arr[0](i, j, k)) * dxInv * mfsq
577  + (flx_v_arr[1](i , j+1, k ) - flx_v_arr[1](i, j, k)) * dyInv * mfsq
578  + (flx_v_arr[2](i , j , k+1) - flx_v_arr[2](i, j, k)) * dzInv );
579  });
580 
581  } else if (v_type == FabType::singlevalued) {
582 
583  ParallelFor(bxy, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
584 
585  if (v_vfrac(i,j,k)>zero) {
586  Real mfsq = mf_vx(i,j,0) * mf_vy(i,j,0);
587 
588  if (v_cflag(i,j,k).isCovered()) {
589 
590  rho_v_rhs(i, j, k) = zero;
591 
592  } else if (v_cflag(i,j,k).isRegular()) {
593 
594  rho_v_rhs(i, j, k) = - ( (v_afrac_x(i+1, j , k ) * flx_v_arr[0](i+1, j , k ) - v_afrac_x(i, j, k) * flx_v_arr[0](i, j, k)) * dxInv * mfsq
595  + (v_afrac_y(i , j+1, k ) * flx_v_arr[1](i , j+1, k ) - v_afrac_y(i, j, k) * flx_v_arr[1](i, j, k)) * dyInv * mfsq
596  + (v_afrac_z(i , j , k+1) * flx_v_arr[2](i , j , k+1) - v_afrac_z(i, j, k) * flx_v_arr[2](i, j, k)) * dzInv ) / v_vfrac(i,j,k);
597  } else {
598  // Bilinear interpolation
599  Real fxm = flx_v_arr[0](i,j,k);
600  if (v_afrac_x(i,j,k) != zero && v_afrac_x(i,j,k) != one) {
601  int jj = j + static_cast<int>(std::copysign(one, v_fcx(i,j,k,0)));
602  int kk = k + static_cast<int>(std::copysign(one, v_fcx(i,j,k,1)));
603  Real fracy = (v_mask(i-1,jj,k) || v_mask(i,jj,k)) ? std::abs(v_fcx(i,j,k,0)) : zero;
604  Real fracz = (v_mask(i-1,j,kk) || v_mask(i,j,kk)) ? std::abs(v_fcx(i,j,k,1)) : zero;
605  fxm = (one-fracy)*(one-fracz)*fxm
606  + fracy *(one-fracz)*flx_v_arr[0](i,jj,k )
607  + fracz *(one-fracy)*flx_v_arr[0](i,j ,kk)
608  + fracy * fracz *flx_v_arr[0](i,jj,kk);
609  }
610 
611  Real fxp = flx_v_arr[0](i+1,j,k);
612  if (v_afrac_x(i+1,j,k) != zero && v_afrac_x(i+1,j,k) != one) {
613  int jj = j + static_cast<int>(std::copysign(one,v_fcx(i+1,j,k,0)));
614  int kk = k + static_cast<int>(std::copysign(one,v_fcx(i+1,j,k,1)));
615  Real fracy = (v_mask(i,jj,k) || v_mask(i+1,jj,k)) ? std::abs(v_fcx(i+1,j,k,0)) : zero;
616  Real fracz = (v_mask(i,j,kk) || v_mask(i+1,j,kk)) ? std::abs(v_fcx(i+1,j,k,1)) : zero;
617  fxp = (one-fracy)*(one-fracz)*fxp
618  + fracy *(one-fracz)*flx_v_arr[0](i+1,jj,k )
619  + fracz *(one-fracy)*flx_v_arr[0](i+1,j ,kk)
620  + fracy * fracz *flx_v_arr[0](i+1,jj,kk);
621  }
622 
623  Real fym = flx_v_arr[1](i,j,k);
624  if (v_afrac_y(i,j,k) != zero && v_afrac_y(i,j,k) != one) {
625  int ii = i + static_cast<int>(std::copysign(one,v_fcy(i,j,k,0)));
626  int kk = k + static_cast<int>(std::copysign(one,v_fcy(i,j,k,1)));
627  Real fracx = (v_mask(ii,j-1,k) || v_mask(ii,j,k)) ? std::abs(v_fcy(i,j,k,0)) : zero;
628  Real fracz = (v_mask(i,j-1,kk) || v_mask(i,j,kk)) ? std::abs(v_fcy(i,j,k,1)) : zero;
629  fym = (one-fracx)*(one-fracz)*fym
630  + fracx *(one-fracz)*flx_v_arr[1](ii,j,k )
631  + fracz *(one-fracx)*flx_v_arr[1](i ,j,kk)
632  + fracx * fracz *flx_v_arr[1](ii,j,kk);
633  }
634 
635  Real fyp = flx_v_arr[1](i,j+1,k);
636  if (v_afrac_y(i,j+1,k) != zero && v_afrac_y(i,j+1,k) != one) {
637  int ii = i + static_cast<int>(std::copysign(one,v_fcy(i,j+1,k,0)));
638  int kk = k + static_cast<int>(std::copysign(one,v_fcy(i,j+1,k,1)));
639  Real fracx = (v_mask(ii,j,k) || v_mask(ii,j+1,k)) ? std::abs(v_fcy(i,j+1,k,0)) : zero;
640  Real fracz = (v_mask(i,j,kk) || v_mask(i,j+1,kk)) ? std::abs(v_fcy(i,j+1,k,1)) : zero;
641  fyp = (one-fracx)*(one-fracz)*fyp
642  + fracx *(one-fracz)*flx_v_arr[1](ii,j+1,k )
643  + fracz *(one-fracx)*flx_v_arr[1](i ,j+1,kk)
644  + fracx * fracz *flx_v_arr[1](ii,j+1,kk);
645  }
646 
647  Real fzm = flx_v_arr[2](i,j,k);
648  if (v_afrac_z(i,j,k) != zero && v_afrac_z(i,j,k) != one) {
649  int ii = i + static_cast<int>(std::copysign(one,v_fcz(i,j,k,0)));
650  int jj = j + static_cast<int>(std::copysign(one,v_fcz(i,j,k,1)));
651  Real fracx = (v_mask(ii,j,k-1) || v_mask(ii,j,k)) ? std::abs(v_fcz(i,j,k,0)) : zero;
652  Real fracy = (v_mask(i,jj,k-1) || v_mask(i,jj,k)) ? std::abs(v_fcz(i,j,k,1)) : zero;
653  fzm = (one-fracx)*(one-fracy)*fzm
654  + fracx *(one-fracy)*flx_v_arr[2](ii,j ,k)
655  + fracy *(one-fracx)*flx_v_arr[2](i ,jj,k)
656  + fracx * fracy *flx_v_arr[2](ii,jj,k);
657  }
658 
659  Real fzp = flx_v_arr[2](i,j,k+1);
660  if (v_afrac_z(i,j,k+1) != zero && v_afrac_z(i,j,k+1) != one) {
661  int ii = i + static_cast<int>(std::copysign(one,v_fcz(i,j,k+1,0)));
662  int jj = j + static_cast<int>(std::copysign(one,v_fcz(i,j,k+1,1)));
663  Real fracx = (v_mask(ii,j,k) || v_mask(ii,j,k+1)) ? std::abs(v_fcz(i,j,k+1,0)) : zero;
664  Real fracy = (v_mask(i,jj,k) || v_mask(i,jj,k+1)) ? std::abs(v_fcz(i,j,k+1,1)) : zero;
665  fzp = (one-fracx)*(one-fracy)*fzp
666  + fracx *(one-fracy)*flx_v_arr[2](ii,j ,k+1)
667  + fracy *(one-fracx)*flx_v_arr[2](i ,jj,k+1)
668  + fracx * fracy *flx_v_arr[2](ii,jj,k+1);
669  }
670 
671  rho_v_rhs(i, j, k) = - ( (v_afrac_x(i+1, j , k ) * fxp - v_afrac_x(i, j, k) * fxm) * dxInv * mfsq
672  + (v_afrac_y(i , j+1, k ) * fyp - v_afrac_y(i, j, k) * fym) * dyInv * mfsq
673  + (v_afrac_z(i , j , k+1) * fzp - v_afrac_z(i, j, k) * fzm) * dzInv ) / v_vfrac(i,j,k);
674  }
675 
676  } else {
677  rho_v_rhs(i, j, k) = zero;
678  }
679  });
680  } // v_type
681 
682  if (w_type == FabType::covered) {
683 
684  ParallelFor(bxz, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
685  rho_w_rhs(i, j, k) = zero;
686  });
687 
688  } else if (w_type == FabType::regular) {
689 
690  ParallelFor(bxz, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
691  Real mfsq = mf_mx(i,j,0) * mf_my(i,j,0);
692  rho_w_rhs(i, j, k) = - ( (flx_w_arr[0](i+1, j , k ) - flx_w_arr[0](i, j, k)) * dxInv * mfsq
693  + (flx_w_arr[1](i , j+1, k ) - flx_w_arr[1](i, j, k)) * dyInv * mfsq
694  + (flx_w_arr[2](i , j , k+1) - flx_w_arr[2](i, j, k)) * dzInv );
695  });
696 
697  } else if (w_type == FabType::singlevalued) {
698 
699  ParallelFor(bxz, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
700  if (w_vfrac(i,j,k)>zero) {
701  Real mfsq = mf_mx(i,j,0) * mf_my(i,j,0);
702 
703  if (w_cflag(i,j,k).isCovered())
704  {
705  rho_w_rhs(i, j, k) = zero;
706  }
707  else if (w_cflag(i,j,k).isRegular())
708  {
709  rho_w_rhs(i, j, k) = - ( (w_afrac_x(i+1, j , k ) * flx_w_arr[0](i+1, j , k ) - w_afrac_x(i, j, k) * flx_w_arr[0](i, j, k)) * dxInv * mfsq
710  + (w_afrac_y(i , j+1, k ) * flx_w_arr[1](i , j+1, k ) - w_afrac_y(i, j, k) * flx_w_arr[1](i, j, k)) * dyInv * mfsq
711  + (w_afrac_z(i , j , k+1) * flx_w_arr[2](i , j , k+1) - w_afrac_z(i, j, k) * flx_w_arr[2](i, j, k)) * dzInv ) / w_vfrac(i,j,k);
712  }
713  else
714  {
715  // Bilinear interpolation
716  Real fxm = flx_w_arr[0](i,j,k);
717  if (w_afrac_x(i,j,k) != zero && w_afrac_x(i,j,k) != one) {
718  int jj = j + static_cast<int>(std::copysign(one, w_fcx(i,j,k,0)));
719  int kk = k + static_cast<int>(std::copysign(one, w_fcx(i,j,k,1)));
720  Real fracy = (w_mask(i-1,jj,k) || w_mask(i,jj,k)) ? std::abs(w_fcx(i,j,k,0)) : zero;
721  Real fracz = (w_mask(i-1,j,kk) || w_mask(i,j,kk)) ? std::abs(w_fcx(i,j,k,1)) : zero;
722  fxm = (one-fracy)*(one-fracz)*fxm
723  + fracy *(one-fracz)*flx_w_arr[0](i,jj,k )
724  + fracz *(one-fracy)*flx_w_arr[0](i,j ,kk)
725  + fracy * fracz *flx_w_arr[0](i,jj,kk);
726  }
727 
728  Real fxp = flx_w_arr[0](i+1,j,k);
729  if (w_afrac_x(i+1,j,k) != zero && w_afrac_x(i+1,j,k) != one) {
730  int jj = j + static_cast<int>(std::copysign(one,w_fcx(i+1,j,k,0)));
731  int kk = k + static_cast<int>(std::copysign(one,w_fcx(i+1,j,k,1)));
732  Real fracy = (w_mask(i,jj,k) || w_mask(i+1,jj,k)) ? std::abs(w_fcx(i+1,j,k,0)) : zero;
733  Real fracz = (w_mask(i,j,kk) || w_mask(i+1,j,kk)) ? std::abs(w_fcx(i+1,j,k,1)) : zero;
734  fxp = (one-fracy)*(one-fracz)*fxp
735  + fracy *(one-fracz)*flx_w_arr[0](i+1,jj,k )
736  + fracz *(one-fracy)*flx_w_arr[0](i+1,j ,kk)
737  + fracy * fracz *flx_w_arr[0](i+1,jj,kk);
738  }
739 
740  Real fym = flx_w_arr[1](i,j,k);
741  if (w_afrac_y(i,j,k) != zero && w_afrac_y(i,j,k) != one) {
742  int ii = i + static_cast<int>(std::copysign(one,w_fcy(i,j,k,0)));
743  int kk = k + static_cast<int>(std::copysign(one,w_fcy(i,j,k,1)));
744  Real fracx = (w_mask(ii,j-1,k) || w_mask(ii,j,k)) ? std::abs(w_fcy(i,j,k,0)) : zero;
745  Real fracz = (w_mask(i,j-1,kk) || w_mask(i,j,kk)) ? std::abs(w_fcy(i,j,k,1)) : zero;
746  fym = (one-fracx)*(one-fracz)*fym
747  + fracx *(one-fracz)*flx_w_arr[1](ii,j,k )
748  + fracz *(one-fracx)*flx_w_arr[1](i ,j,kk)
749  + fracx * fracz *flx_w_arr[1](ii,j,kk);
750  }
751 
752  Real fyp = flx_w_arr[1](i,j+1,k);
753  if (w_afrac_y(i,j+1,k) != zero && w_afrac_y(i,j+1,k) != one) {
754  int ii = i + static_cast<int>(std::copysign(one,w_fcy(i,j+1,k,0)));
755  int kk = k + static_cast<int>(std::copysign(one,w_fcy(i,j+1,k,1)));
756  Real fracx = (w_mask(ii,j,k) || w_mask(ii,j+1,k)) ? std::abs(w_fcy(i,j+1,k,0)) : zero;
757  Real fracz = (w_mask(i,j,kk) || w_mask(i,j+1,kk)) ? std::abs(w_fcy(i,j+1,k,1)) : zero;
758  fyp = (one-fracx)*(one-fracz)*fyp
759  + fracx *(one-fracz)*flx_w_arr[1](ii,j+1,k )
760  + fracz *(one-fracx)*flx_w_arr[1](i ,j+1,kk)
761  + fracx * fracz *flx_w_arr[1](ii,j+1,kk);
762  }
763 
764  Real fzm = flx_w_arr[2](i,j,k);
765  if (w_afrac_z(i,j,k) != zero && w_afrac_z(i,j,k) != one) {
766  int ii = i + static_cast<int>(std::copysign(one,w_fcz(i,j,k,0)));
767  int jj = j + static_cast<int>(std::copysign(one,w_fcz(i,j,k,1)));
768  Real fracx = (w_mask(ii,j,k-1) || w_mask(ii,j,k)) ? std::abs(w_fcz(i,j,k,0)) : zero;
769  Real fracy = (w_mask(i,jj,k-1) || w_mask(i,jj,k)) ? std::abs(w_fcz(i,j,k,1)) : zero;
770  fzm = (one-fracx)*(one-fracy)*fzm
771  + fracx *(one-fracy)*flx_w_arr[2](ii,j ,k)
772  + fracy *(one-fracx)*flx_w_arr[2](i ,jj,k)
773  + fracx * fracy *flx_w_arr[2](ii,jj,k);
774  }
775 
776  Real fzp = flx_w_arr[2](i,j,k+1);
777  if (w_afrac_z(i,j,k+1) != zero && w_afrac_z(i,j,k+1) != one) {
778  int ii = i + static_cast<int>(std::copysign(one,w_fcz(i,j,k+1,0)));
779  int jj = j + static_cast<int>(std::copysign(one,w_fcz(i,j,k+1,1)));
780  Real fracx = (w_mask(ii,j,k) || w_mask(ii,j,k+1)) ? std::abs(w_fcz(i,j,k+1,0)) : zero;
781  Real fracy = (w_mask(i,jj,k) || w_mask(i,jj,k+1)) ? std::abs(w_fcz(i,j,k+1,1)) : zero;
782  fzp = (one-fracx)*(one-fracy)*fzp
783  + fracx *(one-fracy)*flx_w_arr[2](ii,j ,k+1)
784  + fracy *(one-fracx)*flx_w_arr[2](i ,jj,k+1)
785  + fracx * fracy *flx_w_arr[2](ii,jj,k+1);
786  }
787 
788  rho_w_rhs(i, j, k) = - ( (w_afrac_x(i+1, j , k ) * fxp - w_afrac_x(i, j, k) * fxm) * dxInv * mfsq
789  + (w_afrac_y(i , j+1, k ) * fyp - w_afrac_y(i, j, k) * fym) * dyInv * mfsq
790  + (w_afrac_z(i , j , k+1) * fzp - w_afrac_z(i, j, k) * fzm) * dzInv ) / w_vfrac(i,j,k);
791  }
792 
793  } else {
794  rho_w_rhs(i, j, k) = 0;
795  }
796  });
797  } // w_type
798 
799 
800  }
801 
802 }
AdvectionSrcForMom_EB
void AdvectionSrcForMom_EB(const MFIter &mfi, const Box &bxx, const Box &bxy, const Box &bxz, const Vector< Box > &bxx_grown, const Vector< Box > &bxy_grown, const Vector< Box > &bxz_grown, const Array4< Real > &rho_u_rhs, const Array4< Real > &rho_v_rhs, const Array4< Real > &rho_w_rhs, const Array4< const Real > &u, const Array4< const Real > &v, const Array4< const Real > &w, const Array4< const Real > &rho_u, const Array4< const Real > &rho_v, const Array4< const Real > &omega, const GpuArray< Real, AMREX_SPACEDIM > &cellSizeInv, 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 AdvType horiz_adv_type, const AdvType vert_adv_type, const Real horiz_upw_frac, const Real vert_upw_frac, const eb_ &ebfact, GpuArray< Array4< Real >, AMREX_SPACEDIM > &flx_u_arr, GpuArray< Array4< Real >, AMREX_SPACEDIM > &flx_v_arr, GpuArray< Array4< Real >, AMREX_SPACEDIM > &flx_w_arr, const Vector< iMultiFab > &physbnd_mask, const bool already_on_centroids, const int lo_z_face, const int hi_z_face, const Box &)
Definition: ERF_AdvectionSrcForMom_EB.cpp:51
AdvType::Upwind_3rd
@ Upwind_3rd
AdvType::Centered_4th
@ Centered_4th
AdvType::Upwind_5th
@ Upwind_5th
AdvType::Centered_6th
@ Centered_6th
AdvType::Centered_2nd
@ Centered_2nd
dxInv
amrex::GpuArray< Real, AMREX_SPACEDIM > dxInv
Definition: ERF_InitCustomPertVels_ParticleTests.H:17
AMREX_ALWAYS_ASSERT
AMREX_ALWAYS_ASSERT(bx.length()[2]==khi+1)
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);})
one
constexpr amrex::Real one
Definition: ERF_NumericalConstants.H:37
fourth
constexpr amrex::Real fourth
Definition: ERF_NumericalConstants.H:42
zero
constexpr amrex::Real zero
Definition: ERF_NumericalConstants.H:36
w
Real w
Definition: ERF_Plotfile2DInterpolator.cpp:23
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:19
AMREX_ASSERT_WITH_MESSAGE
AMREX_ASSERT_WITH_MESSAGE(wbar_cutoff_min > wbar_cutoff_max, "ERROR: wbar_cutoff_min < wbar_cutoff_max")
eb_::get_w_const_factory
amrex::EBFArrayBoxFactory const * get_w_const_factory() const noexcept
Return the native AMReX z-face EB factory.
Definition: ERF_EB.H:148
eb_::get_u_const_factory
amrex::EBFArrayBoxFactory const * get_u_const_factory() const noexcept
Return the native AMReX x-face EB factory.
Definition: ERF_EB.H:144
eb_::get_v_const_factory
amrex::EBFArrayBoxFactory const * get_v_const_factory() const noexcept
Return the native AMReX y-face EB factory.
Definition: ERF_EB.H:146
IntVars::ymom
@ ymom
Definition: ERF_IndexDefines.H:274
IntVars::zmom
@ zmom
Definition: ERF_IndexDefines.H:275
IntVars::xmom
@ xmom
Definition: ERF_IndexDefines.H:273
MicVar_Morr::omega
@ omega
Definition: ERF_Morrison.H:57
Here is the call graph for this function: