ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
Functions
ERF_EBAdvectionSrcForState.cpp File Reference
#include <ERF_EBAdvection.H>
#include <ERF_EBAdvectionSrcForScalars.H>
#include <ERF_IndexDefines.H>
#include <ERF_TerrainMetrics.H>
#include <ERF_Advection.H>
#include <ERF_AdvectionSrcForScalars.H>
Include dependency graph for ERF_EBAdvectionSrcForState.cpp:

Functions

void EBAdvectionSrcForRho (const Box &bx, const Array4< Real > &advectionSrc, const Array4< const Real > &rho_u, const Array4< const Real > &rho_v, const Array4< const Real > &Omega, const Array4< Real > &avg_xmom, const Array4< Real > &avg_ymom, const Array4< Real > &avg_zmom, const Array4< const int > &mask_arr, const Array4< const EBCellFlag > &cfg_arr, const Array4< const Real > &ax_arr, const Array4< const Real > &ay_arr, const Array4< const Real > &az_arr, const Array4< const Real > &fcx_arr, const Array4< const Real > &fcy_arr, const Array4< const Real > &fcz_arr, const Array4< const Real > &detJ, const GpuArray< Real, AMREX_SPACEDIM > &cellSizeInv, const Array4< const Real > &mf_mx, const Array4< const Real > &mf_my, const Array4< const Real > &mf_uy, const Array4< const Real > &mf_vx, const GpuArray< const Array4< Real >, AMREX_SPACEDIM > &flx_arr, const bool fixed_rho, bool already_on_centroids)
 
void EBAdvectionSrcForScalars (const Box &bx, const int icomp, const int ncomp, const Array4< const Real > &avg_xmom, const Array4< const Real > &avg_ymom, const Array4< const Real > &avg_zmom, const Array4< const Real > &cell_prim, const Array4< Real > &advectionSrc, const Array4< const int > &mask_arr, const Array4< const EBCellFlag > &cfg_arr, const Array4< const Real > &ax_arr, const Array4< const Real > &ay_arr, const Array4< const Real > &az_arr, const Array4< const Real > &fcx_arr, const Array4< const Real > &fcy_arr, const Array4< const Real > &fcz_arr, const Array4< const Real > &detJ, const GpuArray< Real, AMREX_SPACEDIM > &cellSizeInv, const Array4< const Real > &mf_mx, const Array4< const Real > &mf_my, const AdvType horiz_adv_type, const AdvType vert_adv_type, const Real horiz_upw_frac, const Real vert_upw_frac, const GpuArray< const Array4< Real >, AMREX_SPACEDIM > &flx_arr, const Box &domain, const BCRec *bc_ptr_h, bool already_on_centroids)
 

Function Documentation

◆ EBAdvectionSrcForRho()

void EBAdvectionSrcForRho ( const Box &  bx,
const Array4< Real > &  advectionSrc,
const Array4< const Real > &  rho_u,
const Array4< const Real > &  rho_v,
const Array4< const Real > &  Omega,
const Array4< Real > &  avg_xmom,
const Array4< Real > &  avg_ymom,
const Array4< Real > &  avg_zmom,
const Array4< const int > &  mask_arr,
const Array4< const EBCellFlag > &  cfg_arr,
const Array4< const Real > &  ax_arr,
const Array4< const Real > &  ay_arr,
const Array4< const Real > &  az_arr,
const Array4< const Real > &  fcx_arr,
const Array4< const Real > &  fcy_arr,
const Array4< const Real > &  fcz_arr,
const Array4< const Real > &  detJ,
const GpuArray< Real, AMREX_SPACEDIM > &  cellSizeInv,
const Array4< const Real > &  mf_mx,
const Array4< const Real > &  mf_my,
const Array4< const Real > &  mf_uy,
const Array4< const Real > &  mf_vx,
const GpuArray< const Array4< Real >, AMREX_SPACEDIM > &  flx_arr,
const bool  fixed_rho,
bool  already_on_centroids 
)

Function for computing the advective tendency for the update equations for rho and (rho theta) This routine has explicit expressions for all cases (terrain or not) when the horizontal and vertical spatial orders are <= 2, and calls more specialized functions when either (or both) spatial order(s) is greater than 2.

Parameters
[in]bxbox over which the scalars are updated
[out]advectionSrctendency for the scalar update equation
[in]rho_ux-component of momentum
[in]rho_vy-component of momentum
[in]Omegacomponent of momentum normal to the z-coordinate surface
[out]avg_xmomx-component of time-averaged momentum defined in this routine
[out]avg_ymomy-component of time-averaged momentum defined in this routine
[out]avg_zmomz-component of time-averaged momentum defined in this routine
[in]detJJacobian of the metric transformation (= 1 if use_terrain is false)
[in]cellSizeInvinverse of the mesh spacing
[in]mf_mmap factor at cell centers
[in]mf_umap factor at x-faces
[in]mf_vmap factor at y-faces
57 {
58  BL_PROFILE_VAR("EBAdvectionSrcForRho", EBAdvectionSrcForRho);
59  auto dxInv = cellSizeInv[0], dyInv = cellSizeInv[1], dzInv = cellSizeInv[2];
60 
61  const Box xbx = surroundingNodes(bx,0).grow(IntVect(0, 1, 1));
62  const Box ybx = surroundingNodes(bx,1).grow(IntVect(1, 0, 1));
63  const Box zbx = surroundingNodes(bx,2).grow(IntVect(1, 1, 0));
64 
65  ParallelFor(xbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
66  {
67  flx_arr[0](i,j,k,0) = rho_u(i,j,k) / mf_uy(i,j,0);
68  avg_xmom(i,j,k) = flx_arr[0](i,j,k,0);
69  });
70  ParallelFor(ybx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
71  {
72  flx_arr[1](i,j,k,0) = rho_v(i,j,k) / mf_vx(i,j,0);
73  avg_ymom(i,j,k) = flx_arr[1](i,j,k,0);
74  });
75  ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
76  {
77  Real mfsq = mf_mx(i,j,0) * mf_my(i,j,0);
78  flx_arr[2](i,j,k,0) = Omega(i,j,k) / mfsq;
79  avg_zmom(i,j,k) = flx_arr[2](i,j,k,0);
80  });
81 
82  if (fixed_rho) {
83  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
84  {
85  advectionSrc(i,j,k,0) = 0.0;
86  });
87  } else
88  {
89  if (already_on_centroids) {
90 
91  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
92  {
93  if (detJ(i,j,k) > 0.) {
94  Real mfsq = mf_mx(i,j,0) * mf_my(i,j,0);
95  advectionSrc(i,j,k,0) = - mfsq / detJ(i,j,k) * (
96  ( ax_arr(i+1,j,k) * flx_arr[0](i+1,j,k,0) - ax_arr(i,j,k) * flx_arr[0](i,j,k,0) ) * dxInv +
97  ( ay_arr(i,j+1,k) * flx_arr[1](i,j+1,k,0) - ay_arr(i,j,k) * flx_arr[1](i,j,k,0) ) * dyInv +
98  ( az_arr(i,j,k+1) * flx_arr[2](i,j,k+1,0) - az_arr(i,j,k) * flx_arr[2](i,j,k,0) ) * dzInv );
99  } else {
100  advectionSrc(i,j,k,0) = 0.;
101  }
102  });
103 
104  } else {
105 
106  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
107  {
108  if (detJ(i,j,k) > 0.) {
109  Real mfsq = mf_mx(i,j,0) * mf_my(i,j,0);
110  if (cfg_arr(i,j,k).isCovered())
111  {
112  advectionSrc(i,j,k,0) = 0.;
113  }
114  else if (cfg_arr(i,j,k).isRegular())
115  {
116  advectionSrc(i,j,k,0) = - mfsq / detJ(i,j,k) * (
117  ( ax_arr(i+1,j,k) * flx_arr[0](i+1,j,k,0) - ax_arr(i,j,k) * flx_arr[0](i,j,k,0) ) * dxInv +
118  ( ay_arr(i,j+1,k) * flx_arr[1](i,j+1,k,0) - ay_arr(i,j,k) * flx_arr[1](i,j,k,0) ) * dyInv +
119  ( az_arr(i,j,k+1) * flx_arr[2](i,j,k+1,0) - az_arr(i,j,k) * flx_arr[2](i,j,k,0) ) * dzInv );
120  }
121  else
122  {
123  // Bilinear interpolation
124  Real fxm = flx_arr[0](i,j,k,0);
125  if (ax_arr(i,j,k) != Real(0.0) && ax_arr(i,j,k) != Real(1.0)) {
126  int jj = j + static_cast<int>(std::copysign(Real(1.0), fcx_arr(i,j,k,0)));
127  int kk = k + static_cast<int>(std::copysign(Real(1.0), fcx_arr(i,j,k,1)));
128  Real fracy = (mask_arr(i-1,jj,k) || mask_arr(i,jj,k)) ? std::abs(fcx_arr(i,j,k,0)) : Real(0.0);
129  Real fracz = (mask_arr(i-1,j,kk) || mask_arr(i,j,kk)) ? std::abs(fcx_arr(i,j,k,1)) : Real(0.0);
130  fxm = (Real(1.0)-fracy)*(Real(1.0)-fracz)*fxm
131  + fracy *(Real(1.0)-fracz)*flx_arr[0](i,jj,k ,0)
132  + fracz *(Real(1.0)-fracy)*flx_arr[0](i,j ,kk,0)
133  + fracy * fracz *flx_arr[0](i,jj,kk,0);
134  }
135 
136  Real fxp = flx_arr[0](i+1,j,k,0);
137  if (ax_arr(i+1,j,k) != Real(0.0) && ax_arr(i+1,j,k) != Real(1.0)) {
138  int jj = j + static_cast<int>(std::copysign(Real(1.0),fcx_arr(i+1,j,k,0)));
139  int kk = k + static_cast<int>(std::copysign(Real(1.0),fcx_arr(i+1,j,k,1)));
140  Real fracy = (mask_arr(i,jj,k) || mask_arr(i+1,jj,k)) ? std::abs(fcx_arr(i+1,j,k,0)) : Real(0.0);
141  Real fracz = (mask_arr(i,j,kk) || mask_arr(i+1,j,kk)) ? std::abs(fcx_arr(i+1,j,k,1)) : Real(0.0);
142  fxp = (Real(1.0)-fracy)*(Real(1.0)-fracz)*fxp
143  + fracy *(Real(1.0)-fracz)*flx_arr[0](i+1,jj,k ,0)
144  + fracz *(Real(1.0)-fracy)*flx_arr[0](i+1,j ,kk,0)
145  + fracy * fracz *flx_arr[0](i+1,jj,kk,0);
146  }
147 
148  Real fym = flx_arr[1](i,j,k,0);
149  if (ay_arr(i,j,k) != Real(0.0) && ay_arr(i,j,k) != Real(1.0)) {
150  int ii = i + static_cast<int>(std::copysign(Real(1.0),fcy_arr(i,j,k,0)));
151  int kk = k + static_cast<int>(std::copysign(Real(1.0),fcy_arr(i,j,k,1)));
152  Real fracx = (mask_arr(ii,j-1,k) || mask_arr(ii,j,k)) ? std::abs(fcy_arr(i,j,k,0)) : Real(0.0);
153  Real fracz = (mask_arr(i,j-1,kk) || mask_arr(i,j,kk)) ? std::abs(fcy_arr(i,j,k,1)) : Real(0.0);
154  fym = (Real(1.0)-fracx)*(Real(1.0)-fracz)*fym
155  + fracx *(Real(1.0)-fracz)*flx_arr[1](ii,j,k ,0)
156  + fracz *(Real(1.0)-fracx)*flx_arr[1](i ,j,kk,0)
157  + fracx * fracz *flx_arr[1](ii,j,kk,0);
158  }
159 
160  Real fyp = flx_arr[1](i,j+1,k,0);
161  if (ay_arr(i,j+1,k) != Real(0.0) && ay_arr(i,j+1,k) != Real(1.0)) {
162  int ii = i + static_cast<int>(std::copysign(Real(1.0),fcy_arr(i,j+1,k,0)));
163  int kk = k + static_cast<int>(std::copysign(Real(1.0),fcy_arr(i,j+1,k,1)));
164  Real fracx = (mask_arr(ii,j,k) || mask_arr(ii,j+1,k)) ? std::abs(fcy_arr(i,j+1,k,0)) : Real(0.0);
165  Real fracz = (mask_arr(i,j,kk) || mask_arr(i,j+1,kk)) ? std::abs(fcy_arr(i,j+1,k,1)) : Real(0.0);
166  fyp = (Real(1.0)-fracx)*(Real(1.0)-fracz)*fyp
167  + fracx *(Real(1.0)-fracz)*flx_arr[1](ii,j+1,k ,0)
168  + fracz *(Real(1.0)-fracx)*flx_arr[1](i ,j+1,kk,0)
169  + fracx * fracz *flx_arr[1](ii,j+1,kk,0);
170  }
171 
172  Real fzm = flx_arr[2](i,j,k,0);
173  if (az_arr(i,j,k) != Real(0.0) && az_arr(i,j,k) != Real(1.0)) {
174  int ii = i + static_cast<int>(std::copysign(Real(1.0),fcz_arr(i,j,k,0)));
175  int jj = j + static_cast<int>(std::copysign(Real(1.0),fcz_arr(i,j,k,1)));
176  Real fracx = (mask_arr(ii,j,k-1) || mask_arr(ii,j,k)) ? std::abs(fcz_arr(i,j,k,0)) : Real(0.0);
177  Real fracy = (mask_arr(i,jj,k-1) || mask_arr(i,jj,k)) ? std::abs(fcz_arr(i,j,k,1)) : Real(0.0);
178  fzm = (Real(1.0)-fracx)*(Real(1.0)-fracy)*fzm
179  + fracx *(Real(1.0)-fracy)*flx_arr[2](ii,j ,k,0)
180  + fracy *(Real(1.0)-fracx)*flx_arr[2](i ,jj,k,0)
181  + fracx * fracy *flx_arr[2](ii,jj,k,0);
182  }
183 
184  Real fzp = flx_arr[2](i,j,k+1,0);
185  if (az_arr(i,j,k+1) != Real(0.0) && az_arr(i,j,k+1) != Real(1.0)) {
186  int ii = i + static_cast<int>(std::copysign(Real(1.0),fcz_arr(i,j,k+1,0)));
187  int jj = j + static_cast<int>(std::copysign(Real(1.0),fcz_arr(i,j,k+1,1)));
188  Real fracx = (mask_arr(ii,j,k) || mask_arr(ii,j,k+1)) ? std::abs(fcz_arr(i,j,k+1,0)) : Real(0.0);
189  Real fracy = (mask_arr(i,jj,k) || mask_arr(i,jj,k+1)) ? std::abs(fcz_arr(i,j,k+1,1)) : Real(0.0);
190  fzp = (Real(1.0)-fracx)*(Real(1.0)-fracy)*fzp
191  + fracx *(Real(1.0)-fracy)*flx_arr[2](ii,j ,k+1,0)
192  + fracy *(Real(1.0)-fracx)*flx_arr[2](i ,jj,k+1,0)
193  + fracx * fracy *flx_arr[2](ii,jj,k+1,0);
194  }
195 
196  advectionSrc(i,j,k,0) = - mfsq / detJ(i,j,k) * (
197  ( ax_arr(i+1,j,k) * fxp - ax_arr(i,j,k) * fxm ) * dxInv +
198  ( ay_arr(i,j+1,k) * fyp - ay_arr(i,j,k) * fym ) * dyInv +
199  ( az_arr(i,j,k+1) * fzp - az_arr(i,j,k) * fzm ) * dzInv );
200  }
201  } else {
202  advectionSrc(i,j,k,0) = 0.;
203  }
204  });
205 
206  } // already_on_centroids
207  } // fixed_rho
208 }
EBAdvectionSrcForRho
void EBAdvectionSrcForRho(const Box &bx, const Array4< Real > &advectionSrc, const Array4< const Real > &rho_u, const Array4< const Real > &rho_v, const Array4< const Real > &Omega, const Array4< Real > &avg_xmom, const Array4< Real > &avg_ymom, const Array4< Real > &avg_zmom, const Array4< const int > &mask_arr, const Array4< const EBCellFlag > &cfg_arr, const Array4< const Real > &ax_arr, const Array4< const Real > &ay_arr, const Array4< const Real > &az_arr, const Array4< const Real > &fcx_arr, const Array4< const Real > &fcy_arr, const Array4< const Real > &fcz_arr, const Array4< const Real > &detJ, const GpuArray< Real, AMREX_SPACEDIM > &cellSizeInv, const Array4< const Real > &mf_mx, const Array4< const Real > &mf_my, const Array4< const Real > &mf_uy, const Array4< const Real > &mf_vx, const GpuArray< const Array4< Real >, AMREX_SPACEDIM > &flx_arr, const bool fixed_rho, bool already_on_centroids)
Definition: ERF_EBAdvectionSrcForState.cpp:32

◆ EBAdvectionSrcForScalars()

void EBAdvectionSrcForScalars ( const Box &  bx,
const int  icomp,
const int  ncomp,
const Array4< const Real > &  avg_xmom,
const Array4< const Real > &  avg_ymom,
const Array4< const Real > &  avg_zmom,
const Array4< const Real > &  cell_prim,
const Array4< Real > &  advectionSrc,
const Array4< const int > &  mask_arr,
const Array4< const EBCellFlag > &  cfg_arr,
const Array4< const Real > &  ax_arr,
const Array4< const Real > &  ay_arr,
const Array4< const Real > &  az_arr,
const Array4< const Real > &  fcx_arr,
const Array4< const Real > &  fcy_arr,
const Array4< const Real > &  fcz_arr,
const Array4< const Real > &  detJ,
const GpuArray< Real, AMREX_SPACEDIM > &  cellSizeInv,
const Array4< const Real > &  mf_mx,
const Array4< const Real > &  mf_my,
const AdvType  horiz_adv_type,
const AdvType  vert_adv_type,
const Real  horiz_upw_frac,
const Real  vert_upw_frac,
const GpuArray< const Array4< Real >, AMREX_SPACEDIM > &  flx_arr,
const Box &  domain,
const BCRec *  bc_ptr_h,
bool  already_on_centroids 
)

Function for computing the advective tendency for scalars when terrain_type is EB.

Parameters
[in]bxbox over which the scalars are updated if no external boundary conditions
[in]icompcomponent of first scalar to be updated
[in]ncompnumber of components to be updated
[in]avg_xmomx-component of time-averaged momentum defined in this routine
[in]avg_ymomy-component of time-averaged momentum defined in this routine
[in]avg_zmomz-component of time-averaged momentum defined in this routine
[in]cell_primprimitive form of scalar variables, here only potential temperature theta
[out]advectionSrctendency for the scalar update equation
[in]mask_arrCell-centered masks (=1 otherwise, =0 if physbnd)
[in]cfg_arrCell-centered flags
[in]ax_arrArea fraction of x-face
[in]ay_arrArea fraction of y-face
[in]az_arrArea fraction of z-face
[in]fcx_arrFace centroid of x-face
[in]fcy_arrFace centroid of y-face
[in]fcz_arrFace centroid of z-face
[in]detJJacobian of the metric transformation (= 1 if use_terrain is false)
[in]cellSizeInvinverse of the mesh spacing
[in]mf_mxmap factor at cell centers
[in]mf_mymap factor at cell centers
[in]horiz_adv_typeadvection scheme to be used in horiz. directions for dry scalars
[in]vert_adv_typeadvection scheme to be used in vert. directions for dry scalars
[in]horiz_upw_fracupwinding fraction to be used in horiz. directions for dry scalars (for Blended schemes only)
[in]vert_upw_fracupwinding fraction to be used in vert. directions for dry scalars (for Blended schemes only)
268 {
269  BL_PROFILE_VAR("EBAdvectionSrcForScalars", EBAdvectionSrcForScalars);
270  auto dxInv = cellSizeInv[0], dyInv = cellSizeInv[1], dzInv = cellSizeInv[2];
271 
272  const Box xbx = surroundingNodes(bx,0);
273  const Box ybx = surroundingNodes(bx,1);
274  const Box zbx = surroundingNodes(bx,2);
275 
276  // Open bc will be imposed upon all vars (we only access cons here for simplicity)
277  const bool xlo_open = (bc_ptr_h[BCVars::cons_bc].lo(0) == ERFBCType::open);
278  const bool xhi_open = (bc_ptr_h[BCVars::cons_bc].hi(0) == ERFBCType::open);
279  const bool ylo_open = (bc_ptr_h[BCVars::cons_bc].lo(1) == ERFBCType::open);
280  const bool yhi_open = (bc_ptr_h[BCVars::cons_bc].hi(1) == ERFBCType::open);
281 
282  // Only advection operations in bndry normal direction with OPEN BC
283  Box bx_xlo, bx_xhi, bx_ylo, bx_yhi;
284  if (xlo_open) {
285  if ( bx.smallEnd(0) == domain.smallEnd(0)) { bx_xlo = makeSlab( bx,0,domain.smallEnd(0));}
286  }
287  if (xhi_open) {
288  if ( bx.bigEnd(0) == domain.bigEnd(0)) { bx_xhi = makeSlab( bx,0,domain.bigEnd(0) );}
289  }
290  if (ylo_open) {
291  if ( bx.smallEnd(1) == domain.smallEnd(1)) { bx_ylo = makeSlab( bx,1,domain.smallEnd(1));}
292  }
293  if (yhi_open) {
294  if ( bx.bigEnd(1) == domain.bigEnd(1)) { bx_yhi = makeSlab( bx,1,domain.bigEnd(1) );}
295  }
296 
297  // Inline with 2nd order for efficiency
298  // NOTE: For EB, avg_xmom, avg_ymom, avg_zmom were are weighted by area fractions in AdvectionSrcForRho
299  // The flux is weighted by area fraction after its interpolation.
300  if (horiz_adv_type == AdvType::Centered_2nd && vert_adv_type == AdvType::Centered_2nd)
301  {
302  ParallelFor(xbx, ncomp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
303  {
304  const int cons_index = icomp + n;
305  const int prim_index = cons_index - 1;
306  const Real prim_on_face = 0.5 * (cell_prim(i,j,k,prim_index) + cell_prim(i-1,j,k,prim_index));
307  flx_arr[0](i,j,k,cons_index) = avg_xmom(i,j,k) * prim_on_face;
308  });
309  ParallelFor(ybx, ncomp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
310  {
311  const int cons_index = icomp + n;
312  const int prim_index = cons_index - 1;
313  const Real prim_on_face = 0.5 * (cell_prim(i,j,k,prim_index) + cell_prim(i,j-1,k,prim_index));
314  flx_arr[1](i,j,k,cons_index) = avg_ymom(i,j,k) * prim_on_face;
315  });
316  ParallelFor(zbx, ncomp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
317  {
318  const int cons_index = icomp + n;
319  const int prim_index = cons_index - 1;
320  const Real prim_on_face = 0.5 * (cell_prim(i,j,k,prim_index) + cell_prim(i,j,k-1,prim_index));
321  flx_arr[2](i,j,k,cons_index) = avg_zmom(i,j,k) * prim_on_face;
322  });
323 
324  // Template higher order methods (horizontal first)
325  } else {
326  switch(horiz_adv_type) {
327  case AdvType::Centered_2nd:
328  EBAdvectionSrcForScalarsVert<CENTERED2>(bx, ncomp, icomp, flx_arr, cell_prim,
329  avg_xmom, avg_ymom, avg_zmom, cfg_arr, ax_arr, ay_arr, az_arr,
330  horiz_upw_frac, vert_upw_frac, vert_adv_type);
331  break;
332  case AdvType::Upwind_3rd:
333  EBAdvectionSrcForScalarsVert<UPWIND3>(bx, ncomp, icomp, flx_arr, cell_prim,
334  avg_xmom, avg_ymom, avg_zmom, cfg_arr, ax_arr, ay_arr, az_arr,
335  horiz_upw_frac, vert_upw_frac, vert_adv_type);
336  break;
337  case AdvType::Centered_4th:
338  EBAdvectionSrcForScalarsVert<CENTERED4>(bx, ncomp, icomp, flx_arr, cell_prim,
339  avg_xmom, avg_ymom, avg_zmom, cfg_arr, ax_arr, ay_arr, az_arr,
340  horiz_upw_frac, vert_upw_frac, vert_adv_type);
341  break;
342  case AdvType::Upwind_5th:
343  EBAdvectionSrcForScalarsVert<UPWIND5>(bx, ncomp, icomp, flx_arr, cell_prim,
344  avg_xmom, avg_ymom, avg_zmom, cfg_arr, ax_arr, ay_arr, az_arr,
345  horiz_upw_frac, vert_upw_frac, vert_adv_type);
346  break;
347  case AdvType::Centered_6th:
348  EBAdvectionSrcForScalarsVert<CENTERED6>(bx, ncomp, icomp, flx_arr, cell_prim,
349  avg_xmom, avg_ymom, avg_zmom, cfg_arr, ax_arr, ay_arr, az_arr,
350  horiz_upw_frac, vert_upw_frac, vert_adv_type);
351  break;
352  default:
353  AMREX_ASSERT_WITH_MESSAGE(false, "Unknown advection scheme! WENO is currently not supported for EB.");
354  }
355  }
356 
357  // Compute divergence
358 
359  if (already_on_centroids) {
360 
361  ParallelFor(bx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
362  {
363  const int cons_index = icomp + n;
364  if (detJ(i,j,k) > 0.)
365  {
366  Real invdetJ = 1.0 / detJ(i,j,k);
367  Real mfsq = mf_mx(i,j,0) * mf_my(i,j,0);
368 
369  advectionSrc(i,j,k,cons_index) = - invdetJ * mfsq * (
370  ( ax_arr(i+1,j,k) * flx_arr[0](i+1,j,k,cons_index) - ax_arr(i,j,k) * flx_arr[0](i ,j,k,cons_index) ) * dxInv +
371  ( ay_arr(i,j+1,k) * flx_arr[1](i,j+1,k,cons_index) - ay_arr(i,j,k) * flx_arr[1](i,j ,k,cons_index) ) * dyInv +
372  ( az_arr(i,j,k+1) * flx_arr[2](i,j,k+1,cons_index) - az_arr(i,j,k) * flx_arr[2](i,j,k ,cons_index) ) * dzInv );
373  } else {
374  advectionSrc(i,j,k,cons_index) = 0.;
375  }
376  });
377 
378  } else {
379 
380  AMREX_HOST_DEVICE_FOR_4D(bx,ncomp,i,j,k,n,
381  {
382  const int cons_index = icomp + n;
383  if (detJ(i,j,k) > 0.)
384  {
385  Real invdetJ = 1.0 / detJ(i,j,k);
386  Real mfsq = mf_mx(i,j,0) * mf_my(i,j,0);
387  if (cfg_arr(i,j,k).isCovered())
388  {
389  advectionSrc(i,j,k,cons_index) = Real(0.0);
390  }
391  else if (cfg_arr(i,j,k).isRegular())
392  {
393  advectionSrc(i,j,k,cons_index) = - invdetJ * mfsq * (
394  ( ax_arr(i+1,j,k) * flx_arr[0](i+1,j,k,cons_index) - ax_arr(i,j,k) * flx_arr[0](i ,j,k,cons_index) ) * dxInv +
395  ( ay_arr(i,j+1,k) * flx_arr[1](i,j+1,k,cons_index) - ay_arr(i,j,k) * flx_arr[1](i,j ,k,cons_index) ) * dyInv +
396  ( az_arr(i,j,k+1) * flx_arr[2](i,j,k+1,cons_index) - az_arr(i,j,k) * flx_arr[2](i,j,k ,cons_index) ) * dzInv );
397  }
398  else
399  {
400  // Bilinear interpolation
401  Real fxm = flx_arr[0](i,j,k,cons_index);
402  if (ax_arr(i,j,k) != Real(0.0) && ax_arr(i,j,k) != Real(1.0)) {
403  int jj = j + static_cast<int>(std::copysign(Real(1.0), fcx_arr(i,j,k,0)));
404  int kk = k + static_cast<int>(std::copysign(Real(1.0), fcx_arr(i,j,k,1)));
405  Real fracy = (mask_arr(i-1,jj,k) || mask_arr(i,jj,k)) ? std::abs(fcx_arr(i,j,k,0)) : Real(0.0);
406  Real fracz = (mask_arr(i-1,j,kk) || mask_arr(i,j,kk)) ? std::abs(fcx_arr(i,j,k,1)) : Real(0.0);
407  fxm = (Real(1.0)-fracy)*(Real(1.0)-fracz)*fxm
408  + fracy *(Real(1.0)-fracz)*flx_arr[0](i,jj,k ,cons_index)
409  + fracz *(Real(1.0)-fracy)*flx_arr[0](i,j ,kk,cons_index)
410  + fracy * fracz *flx_arr[0](i,jj,kk,cons_index);
411  }
412 
413  Real fxp = flx_arr[0](i+1,j,k,cons_index);
414  if (ax_arr(i+1,j,k) != Real(0.0) && ax_arr(i+1,j,k) != Real(1.0)) {
415  int jj = j + static_cast<int>(std::copysign(Real(1.0),fcx_arr(i+1,j,k,0)));
416  int kk = k + static_cast<int>(std::copysign(Real(1.0),fcx_arr(i+1,j,k,1)));
417  Real fracy = (mask_arr(i,jj,k) || mask_arr(i+1,jj,k)) ? std::abs(fcx_arr(i+1,j,k,0)) : Real(0.0);
418  Real fracz = (mask_arr(i,j,kk) || mask_arr(i+1,j,kk)) ? std::abs(fcx_arr(i+1,j,k,1)) : Real(0.0);
419  fxp = (Real(1.0)-fracy)*(Real(1.0)-fracz)*fxp
420  + fracy *(Real(1.0)-fracz)*flx_arr[0](i+1,jj,k ,cons_index)
421  + fracz *(Real(1.0)-fracy)*flx_arr[0](i+1,j ,kk,cons_index)
422  + fracy * fracz *flx_arr[0](i+1,jj,kk,cons_index);
423  }
424 
425  Real fym = flx_arr[1](i,j,k,cons_index);
426  if (ay_arr(i,j,k) != Real(0.0) && ay_arr(i,j,k) != Real(1.0)) {
427  int ii = i + static_cast<int>(std::copysign(Real(1.0),fcy_arr(i,j,k,0)));
428  int kk = k + static_cast<int>(std::copysign(Real(1.0),fcy_arr(i,j,k,1)));
429  Real fracx = (mask_arr(ii,j-1,k) || mask_arr(ii,j,k)) ? std::abs(fcy_arr(i,j,k,0)) : Real(0.0);
430  Real fracz = (mask_arr(i,j-1,kk) || mask_arr(i,j,kk)) ? std::abs(fcy_arr(i,j,k,1)) : Real(0.0);
431  fym = (Real(1.0)-fracx)*(Real(1.0)-fracz)*fym
432  + fracx *(Real(1.0)-fracz)*flx_arr[1](ii,j,k ,cons_index)
433  + fracz *(Real(1.0)-fracx)*flx_arr[1](i ,j,kk,cons_index)
434  + fracx * fracz *flx_arr[1](ii,j,kk,cons_index);
435  }
436 
437  Real fyp = flx_arr[1](i,j+1,k,cons_index);
438  if (ay_arr(i,j+1,k) != Real(0.0) && ay_arr(i,j+1,k) != Real(1.0)) {
439  int ii = i + static_cast<int>(std::copysign(Real(1.0),fcy_arr(i,j+1,k,0)));
440  int kk = k + static_cast<int>(std::copysign(Real(1.0),fcy_arr(i,j+1,k,1)));
441  Real fracx = (mask_arr(ii,j,k) || mask_arr(ii,j+1,k)) ? std::abs(fcy_arr(i,j+1,k,0)) : Real(0.0);
442  Real fracz = (mask_arr(i,j,kk) || mask_arr(i,j+1,kk)) ? std::abs(fcy_arr(i,j+1,k,1)) : Real(0.0);
443  fyp = (Real(1.0)-fracx)*(Real(1.0)-fracz)*fyp
444  + fracx *(Real(1.0)-fracz)*flx_arr[1](ii,j+1,k ,cons_index)
445  + fracz *(Real(1.0)-fracx)*flx_arr[1](i ,j+1,kk,cons_index)
446  + fracx * fracz *flx_arr[1](ii,j+1,kk,cons_index);
447  }
448 
449  Real fzm = flx_arr[2](i,j,k,cons_index);
450  if (az_arr(i,j,k) != Real(0.0) && az_arr(i,j,k) != Real(1.0)) {
451  int ii = i + static_cast<int>(std::copysign(Real(1.0),fcz_arr(i,j,k,0)));
452  int jj = j + static_cast<int>(std::copysign(Real(1.0),fcz_arr(i,j,k,1)));
453  Real fracx = (mask_arr(ii,j,k-1) || mask_arr(ii,j,k)) ? std::abs(fcz_arr(i,j,k,0)) : Real(0.0);
454  Real fracy = (mask_arr(i,jj,k-1) || mask_arr(i,jj,k)) ? std::abs(fcz_arr(i,j,k,1)) : Real(0.0);
455  fzm = (Real(1.0)-fracx)*(Real(1.0)-fracy)*fzm
456  + fracx *(Real(1.0)-fracy)*flx_arr[2](ii,j ,k,cons_index)
457  + fracy *(Real(1.0)-fracx)*flx_arr[2](i ,jj,k,cons_index)
458  + fracx * fracy *flx_arr[2](ii,jj,k,cons_index);
459  }
460 
461  Real fzp = flx_arr[2](i,j,k+1,cons_index);
462  if (az_arr(i,j,k+1) != Real(0.0) && az_arr(i,j,k+1) != Real(1.0)) {
463  int ii = i + static_cast<int>(std::copysign(Real(1.0),fcz_arr(i,j,k+1,0)));
464  int jj = j + static_cast<int>(std::copysign(Real(1.0),fcz_arr(i,j,k+1,1)));
465  Real fracx = (mask_arr(ii,j,k) || mask_arr(ii,j,k+1)) ? std::abs(fcz_arr(i,j,k+1,0)) : Real(0.0);
466  Real fracy = (mask_arr(i,jj,k) || mask_arr(i,jj,k+1)) ? std::abs(fcz_arr(i,j,k+1,1)) : Real(0.0);
467  fzp = (Real(1.0)-fracx)*(Real(1.0)-fracy)*fzp
468  + fracx *(Real(1.0)-fracy)*flx_arr[2](ii,j ,k+1,cons_index)
469  + fracy *(Real(1.0)-fracx)*flx_arr[2](i ,jj,k+1,cons_index)
470  + fracx * fracy *flx_arr[2](ii,jj,k+1,cons_index);
471  }
472 
473  advectionSrc(i,j,k,cons_index) = - invdetJ * mfsq * (
474  ( ax_arr(i+1,j,k) * fxp - ax_arr(i,j,k) * fxm ) * dxInv
475  + ( ay_arr(i,j+1,k) * fyp - ay_arr(i,j,k) * fym ) * dyInv
476  + ( az_arr(i,j,k+1) * fzp - az_arr(i,j,k) * fzm ) * dzInv );
477  }
478 
479  // eb_compute_divergence(i,j,k,n,advectionSrc,AMREX_D_DECL(flx_arr[0],flx_arr[1],flx_arr[2]),
480  // mask_arr, cfg_arr, detJ, AMREX_D_DECL(ax_arr,ay_arr,az_arr),
481  // AMREX_D_DECL(fcx_arr,fcy_arr,fcz_arr), cellSizeInv, already_on_centroids);
482 
483 
484  } else {
485  advectionSrc(i,j,k,cons_index) = 0.;
486  }
487  });
488 
489  }
490 
491  // Special advection operator for open BC (bndry tangent operations)
492  if (xlo_open) {
493  bool do_lo = true;
494  AdvectionSrcForOpenBC_Tangent_Cons(bx_xlo, 0, icomp, ncomp, advectionSrc, cell_prim,
495  avg_xmom, avg_ymom, avg_zmom,
496  detJ, cellSizeInv, do_lo);
497  }
498  if (xhi_open) {
499  AdvectionSrcForOpenBC_Tangent_Cons(bx_xhi, 0, icomp, ncomp, advectionSrc, cell_prim,
500  avg_xmom, avg_ymom, avg_zmom,
501  detJ, cellSizeInv);
502  }
503  if (ylo_open) {
504  bool do_lo = true;
505  AdvectionSrcForOpenBC_Tangent_Cons(bx_ylo, 1, icomp, ncomp, advectionSrc, cell_prim,
506  avg_xmom, avg_ymom, avg_zmom,
507  detJ, cellSizeInv, do_lo);
508  }
509  if (yhi_open) {
510  AdvectionSrcForOpenBC_Tangent_Cons(bx_yhi, 1, icomp, ncomp, advectionSrc, cell_prim,
511  avg_xmom, avg_ymom, avg_zmom,
512  detJ, cellSizeInv);
513  }
514 }
AdvectionSrcForOpenBC_Tangent_Cons
void AdvectionSrcForOpenBC_Tangent_Cons(const amrex::Box &bx, const int &dir, const int &icomp, const int &ncomp, const amrex::Array4< amrex::Real > &cell_rhs, const amrex::Array4< const amrex::Real > &cell_prim, const amrex::Array4< const amrex::Real > &avg_xmom, const amrex::Array4< const amrex::Real > &avg_ymom, const amrex::Array4< const amrex::Real > &avg_zmom, const amrex::Array4< const amrex::Real > &detJ, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const bool do_lo=false)
EBAdvectionSrcForScalars
void EBAdvectionSrcForScalars(const Box &bx, const int icomp, const int ncomp, const Array4< const Real > &avg_xmom, const Array4< const Real > &avg_ymom, const Array4< const Real > &avg_zmom, const Array4< const Real > &cell_prim, const Array4< Real > &advectionSrc, const Array4< const int > &mask_arr, const Array4< const EBCellFlag > &cfg_arr, const Array4< const Real > &ax_arr, const Array4< const Real > &ay_arr, const Array4< const Real > &az_arr, const Array4< const Real > &fcx_arr, const Array4< const Real > &fcy_arr, const Array4< const Real > &fcz_arr, const Array4< const Real > &detJ, const GpuArray< Real, AMREX_SPACEDIM > &cellSizeInv, const Array4< const Real > &mf_mx, const Array4< const Real > &mf_my, const AdvType horiz_adv_type, const AdvType vert_adv_type, const Real horiz_upw_frac, const Real vert_upw_frac, const GpuArray< const Array4< Real >, AMREX_SPACEDIM > &flx_arr, const Box &domain, const BCRec *bc_ptr_h, bool already_on_centroids)
Definition: ERF_EBAdvectionSrcForState.cpp:240
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
BCVars::cons_bc
@ cons_bc
Definition: ERF_IndexDefines.H:76
ERFBCType::open
@ open
Definition: ERF_IndexDefines.H:215
Here is the call graph for this function: