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