ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
Functions
ERF_EBAdvectionSrcForMom.H File Reference

Defines EB momentum-advection flux builders for higher-order stencils. More...

#include <ERF_IndexDefines.H>
#include <ERF_Interpolation.H>
Include dependency graph for ERF_EBAdvectionSrcForMom.H:
This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Functions

template<typename InterpType_H , typename InterpType_V , typename WallInterpType >
void EBAdvectionSrcForMomWrapper (const amrex::Vector< amrex::Box > &bxx_grown, const amrex::Vector< amrex::Box > &bxy_grown, const amrex::Vector< amrex::Box > &bxz_grown, const amrex::Array4< const amrex::Real > &rho_u, const amrex::Array4< const amrex::Real > &rho_v, const amrex::Array4< const amrex::Real > &rho_w, const amrex::Array4< const amrex::Real > &u, const amrex::Array4< const amrex::Real > &v, const amrex::Array4< const amrex::Real > &w, const amrex::FabType u_type, const amrex::FabType v_type, const amrex::FabType w_type, const amrex::Array4< const amrex::EBCellFlag > &u_cflag, const amrex::Array4< const amrex::Real > &u_afrac_x, const amrex::Array4< const amrex::Real > &u_afrac_y, const amrex::Array4< const amrex::Real > &u_afrac_z, const amrex::Array4< const amrex::EBCellFlag > &v_cflag, const amrex::Array4< const amrex::Real > &v_afrac_x, const amrex::Array4< const amrex::Real > &v_afrac_y, const amrex::Array4< const amrex::Real > &v_afrac_z, const amrex::Array4< const amrex::EBCellFlag > &w_cflag, const amrex::Array4< const amrex::Real > &w_afrac_x, const amrex::Array4< const amrex::Real > &w_afrac_y, const amrex::Array4< const amrex::Real > &w_afrac_z, const amrex::Array4< const amrex::Real > &mf_uy_inv, const amrex::Array4< const amrex::Real > &mf_vx_inv, const amrex::Real upw_frac_h, const amrex::Real upw_frac_v, const AdvType, amrex::GpuArray< amrex::Array4< amrex::Real >, AMREX_SPACEDIM > &flx_u_arr, amrex::GpuArray< amrex::Array4< amrex::Real >, AMREX_SPACEDIM > &flx_v_arr, amrex::GpuArray< amrex::Array4< amrex::Real >, AMREX_SPACEDIM > &flx_w_arr, const int lo_z_face, const int hi_z_face)
 Build EB momentum fluxes with selected horizontal and vertical interpolation. More...
 
template<typename InterpType_H >
void EBAdvectionSrcForMomVert (const amrex::Vector< amrex::Box > &bxx_grown, const amrex::Vector< amrex::Box > &bxy_grown, const amrex::Vector< amrex::Box > &bxz_grown, const amrex::Array4< const amrex::Real > &rho_u, const amrex::Array4< const amrex::Real > &rho_v, const amrex::Array4< const amrex::Real > &rho_w, const amrex::Array4< const amrex::Real > &u, const amrex::Array4< const amrex::Real > &v, const amrex::Array4< const amrex::Real > &w, const amrex::FabType u_type, const amrex::FabType v_type, const amrex::FabType w_type, const amrex::Array4< const amrex::EBCellFlag > &u_cflag, const amrex::Array4< const amrex::Real > &u_afrac_x, const amrex::Array4< const amrex::Real > &u_afrac_y, const amrex::Array4< const amrex::Real > &u_afrac_z, const amrex::Array4< const amrex::EBCellFlag > &v_cflag, const amrex::Array4< const amrex::Real > &v_afrac_x, const amrex::Array4< const amrex::Real > &v_afrac_y, const amrex::Array4< const amrex::Real > &v_afrac_z, const amrex::Array4< const amrex::EBCellFlag > &w_cflag, const amrex::Array4< const amrex::Real > &w_afrac_x, const amrex::Array4< const amrex::Real > &w_afrac_y, const amrex::Array4< const amrex::Real > &w_afrac_z, const amrex::Array4< const amrex::Real > &mf_uy_inv, const amrex::Array4< const amrex::Real > &mf_vx_inv, const amrex::Real upw_frac_h, const amrex::Real upw_frac_v, const AdvType vert_adv_type, amrex::GpuArray< amrex::Array4< amrex::Real >, AMREX_SPACEDIM > &flx_u_arr, amrex::GpuArray< amrex::Array4< amrex::Real >, AMREX_SPACEDIM > &flx_v_arr, amrex::GpuArray< amrex::Array4< amrex::Real >, AMREX_SPACEDIM > &flx_w_arr, const int lo_z_face, const int hi_z_face)
 Dispatch EB momentum flux construction on the vertical advection scheme. More...
 

Detailed Description

Defines EB momentum-advection flux builders for higher-order stencils.

Function Documentation

◆ EBAdvectionSrcForMomVert()

template<typename InterpType_H >
void EBAdvectionSrcForMomVert ( const amrex::Vector< amrex::Box > &  bxx_grown,
const amrex::Vector< amrex::Box > &  bxy_grown,
const amrex::Vector< amrex::Box > &  bxz_grown,
const amrex::Array4< const amrex::Real > &  rho_u,
const amrex::Array4< const amrex::Real > &  rho_v,
const amrex::Array4< const amrex::Real > &  rho_w,
const amrex::Array4< const amrex::Real > &  u,
const amrex::Array4< const amrex::Real > &  v,
const amrex::Array4< const amrex::Real > &  w,
const amrex::FabType  u_type,
const amrex::FabType  v_type,
const amrex::FabType  w_type,
const amrex::Array4< const amrex::EBCellFlag > &  u_cflag,
const amrex::Array4< const amrex::Real > &  u_afrac_x,
const amrex::Array4< const amrex::Real > &  u_afrac_y,
const amrex::Array4< const amrex::Real > &  u_afrac_z,
const amrex::Array4< const amrex::EBCellFlag > &  v_cflag,
const amrex::Array4< const amrex::Real > &  v_afrac_x,
const amrex::Array4< const amrex::Real > &  v_afrac_y,
const amrex::Array4< const amrex::Real > &  v_afrac_z,
const amrex::Array4< const amrex::EBCellFlag > &  w_cflag,
const amrex::Array4< const amrex::Real > &  w_afrac_x,
const amrex::Array4< const amrex::Real > &  w_afrac_y,
const amrex::Array4< const amrex::Real > &  w_afrac_z,
const amrex::Array4< const amrex::Real > &  mf_uy_inv,
const amrex::Array4< const amrex::Real > &  mf_vx_inv,
const amrex::Real  upw_frac_h,
const amrex::Real  upw_frac_v,
const AdvType  vert_adv_type,
amrex::GpuArray< amrex::Array4< amrex::Real >, AMREX_SPACEDIM > &  flx_u_arr,
amrex::GpuArray< amrex::Array4< amrex::Real >, AMREX_SPACEDIM > &  flx_v_arr,
amrex::GpuArray< amrex::Array4< amrex::Real >, AMREX_SPACEDIM > &  flx_w_arr,
const int  lo_z_face,
const int  hi_z_face 
)

Dispatch EB momentum flux construction on the vertical advection scheme.

Template Parameters
InterpType_HHorizontal interpolation functor selected by the caller.
612 {
613  if (vert_adv_type == AdvType::Centered_2nd) {
614  EBAdvectionSrcForMomWrapper<InterpType_H,CENTERED2,UPWINDALL>(
615  bxx_grown, bxy_grown, bxz_grown,
616  rho_u, rho_v, rho_w, u, v, w,
617  u_type, v_type, w_type,
618  u_cflag, u_afrac_x, u_afrac_y, u_afrac_z,
619  v_cflag, v_afrac_x, v_afrac_y, v_afrac_z,
620  w_cflag, w_afrac_x, w_afrac_y, w_afrac_z,
621  mf_uy_inv, mf_vx_inv,
622  upw_frac_h, upw_frac_v, vert_adv_type,
623  flx_u_arr, flx_v_arr, flx_w_arr,
624  lo_z_face, hi_z_face);
625  } else if (vert_adv_type == AdvType::Upwind_3rd) {
626  EBAdvectionSrcForMomWrapper<InterpType_H,UPWIND3,UPWINDALL>(
627  bxx_grown, bxy_grown, bxz_grown,
628  rho_u, rho_v, rho_w, u, v, w,
629  u_type, v_type, w_type,
630  u_cflag, u_afrac_x, u_afrac_y, u_afrac_z,
631  v_cflag, v_afrac_x, v_afrac_y, v_afrac_z,
632  w_cflag, w_afrac_x, w_afrac_y, w_afrac_z,
633  mf_uy_inv, mf_vx_inv,
634  upw_frac_h, upw_frac_v, vert_adv_type,
635  flx_u_arr, flx_v_arr, flx_w_arr,
636  lo_z_face, hi_z_face);
637  } else if (vert_adv_type == AdvType::Centered_4th) {
638  EBAdvectionSrcForMomWrapper<InterpType_H,CENTERED4,UPWINDALL>(
639  bxx_grown, bxy_grown, bxz_grown,
640  rho_u, rho_v, rho_w, u, v, w,
641  u_type, v_type, w_type,
642  u_cflag, u_afrac_x, u_afrac_y, u_afrac_z,
643  v_cflag, v_afrac_x, v_afrac_y, v_afrac_z,
644  w_cflag, w_afrac_x, w_afrac_y, w_afrac_z,
645  mf_uy_inv, mf_vx_inv,
646  upw_frac_h, upw_frac_v, vert_adv_type,
647  flx_u_arr, flx_v_arr, flx_w_arr,
648  lo_z_face, hi_z_face);
649  } else if (vert_adv_type == AdvType::Upwind_5th) {
650  EBAdvectionSrcForMomWrapper<InterpType_H,UPWIND5,UPWINDALL>(
651  bxx_grown, bxy_grown, bxz_grown,
652  rho_u, rho_v, rho_w, u, v, w,
653  u_type, v_type, w_type,
654  u_cflag, u_afrac_x, u_afrac_y, u_afrac_z,
655  v_cflag, v_afrac_x, v_afrac_y, v_afrac_z,
656  w_cflag, w_afrac_x, w_afrac_y, w_afrac_z,
657  mf_uy_inv, mf_vx_inv,
658  upw_frac_h, upw_frac_v, vert_adv_type,
659  flx_u_arr, flx_v_arr, flx_w_arr,
660  lo_z_face, hi_z_face);
661  } else if (vert_adv_type == AdvType::Centered_6th) {
662  EBAdvectionSrcForMomWrapper<InterpType_H,CENTERED6,UPWINDALL>(
663  bxx_grown, bxy_grown, bxz_grown,
664  rho_u, rho_v, rho_w, u, v, w,
665  u_type, v_type, w_type,
666  u_cflag, u_afrac_x, u_afrac_y, u_afrac_z,
667  v_cflag, v_afrac_x, v_afrac_y, v_afrac_z,
668  w_cflag, w_afrac_x, w_afrac_y, w_afrac_z,
669  mf_uy_inv, mf_vx_inv,
670  upw_frac_h, upw_frac_v, vert_adv_type,
671  flx_u_arr, flx_v_arr, flx_w_arr,
672  lo_z_face, hi_z_face);
673  } else {
674  AMREX_ASSERT_WITH_MESSAGE(false, "Unknown advection scheme!");
675  }
676 }
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
w
Real w
Definition: ERF_Plotfile2DInterpolator.cpp:23
AMREX_ASSERT_WITH_MESSAGE
AMREX_ASSERT_WITH_MESSAGE(wbar_cutoff_min > wbar_cutoff_max, "ERROR: wbar_cutoff_min < wbar_cutoff_max")
Here is the call graph for this function:

◆ EBAdvectionSrcForMomWrapper()

template<typename InterpType_H , typename InterpType_V , typename WallInterpType >
void EBAdvectionSrcForMomWrapper ( const amrex::Vector< amrex::Box > &  bxx_grown,
const amrex::Vector< amrex::Box > &  bxy_grown,
const amrex::Vector< amrex::Box > &  bxz_grown,
const amrex::Array4< const amrex::Real > &  rho_u,
const amrex::Array4< const amrex::Real > &  rho_v,
const amrex::Array4< const amrex::Real > &  rho_w,
const amrex::Array4< const amrex::Real > &  u,
const amrex::Array4< const amrex::Real > &  v,
const amrex::Array4< const amrex::Real > &  w,
const amrex::FabType  u_type,
const amrex::FabType  v_type,
const amrex::FabType  w_type,
const amrex::Array4< const amrex::EBCellFlag > &  u_cflag,
const amrex::Array4< const amrex::Real > &  u_afrac_x,
const amrex::Array4< const amrex::Real > &  u_afrac_y,
const amrex::Array4< const amrex::Real > &  u_afrac_z,
const amrex::Array4< const amrex::EBCellFlag > &  v_cflag,
const amrex::Array4< const amrex::Real > &  v_afrac_x,
const amrex::Array4< const amrex::Real > &  v_afrac_y,
const amrex::Array4< const amrex::Real > &  v_afrac_z,
const amrex::Array4< const amrex::EBCellFlag > &  w_cflag,
const amrex::Array4< const amrex::Real > &  w_afrac_x,
const amrex::Array4< const amrex::Real > &  w_afrac_y,
const amrex::Array4< const amrex::Real > &  w_afrac_z,
const amrex::Array4< const amrex::Real > &  mf_uy_inv,
const amrex::Array4< const amrex::Real > &  mf_vx_inv,
const amrex::Real  upw_frac_h,
const amrex::Real  upw_frac_v,
const  AdvType,
amrex::GpuArray< amrex::Array4< amrex::Real >, AMREX_SPACEDIM > &  flx_u_arr,
amrex::GpuArray< amrex::Array4< amrex::Real >, AMREX_SPACEDIM > &  flx_v_arr,
amrex::GpuArray< amrex::Array4< amrex::Real >, AMREX_SPACEDIM > &  flx_w_arr,
const int  lo_z_face,
const int  hi_z_face 
)

Build EB momentum fluxes with selected horizontal and vertical interpolation.

The wrapper handles the x-, y-, and z-momentum grids separately and reduces interpolation order near covered EB cells and vertical domain boundaries. It stores unweighted flux densities on open faces and zero on closed faces; area fractions are applied by the divergence after centroid interpolation.

Template Parameters
InterpType_HHorizontal interpolation functor.
InterpType_VVertical interpolation functor.
WallInterpTypeWall interpolation functor retained for interface symmetry.
55 {
56  // Instantiate the appropriate structs
57  InterpType_H interp_u_horz(u, upw_frac_h); InterpType_V interp_u_vert(u, upw_frac_v); // X-MOM
58  InterpType_H interp_v_horz(v, upw_frac_h); InterpType_V interp_v_vert(v, upw_frac_v); // Y-MOM
59  InterpType_H interp_w_horz(w, upw_frac_h); InterpType_V interp_w_vert(w, upw_frac_v); // Z-MOM
60  // WallInterpType interp_w_wall(w, upw_frac_v); // Z-MOM @ wall
61 
62  const int ncells_u_horz = interp_u_horz.GetUpwindCellNumber();
63  const int ncells_u_vert = interp_u_vert.GetUpwindCellNumber();
64  const int ncells_v_horz = interp_v_horz.GetUpwindCellNumber();
65  const int ncells_v_vert = interp_v_vert.GetUpwindCellNumber();
66  const int ncells_w_horz = interp_w_horz.GetUpwindCellNumber();
67  const int ncells_w_vert = interp_w_vert.GetUpwindCellNumber();
68 
69  // Instantiate lower-order interpolation structs with direction-specific upwinding
70  CENTERED2 interp_u_CEN2(u, 0);
71  CENTERED2 interp_v_CEN2(v, 0);
72  CENTERED2 interp_w_CEN2(w, 0);
73  CENTERED4 interp_w_CEN4(w, 0);
74 
75  UPWIND3 interp_u_UPW3_horz(u, upw_frac_h); UPWIND3 interp_u_UPW3_vert(u, upw_frac_v);
76  UPWIND3 interp_v_UPW3_horz(v, upw_frac_h); UPWIND3 interp_v_UPW3_vert(v, upw_frac_v);
77  UPWIND3 interp_w_UPW3_horz(w, upw_frac_h); UPWIND3 interp_w_UPW3_vert(w, upw_frac_v);
78 
79  // The vertical fallbacks take the vertical upwind fraction
80  UPWIND3 interp_u_UPW3_v(u, upw_frac_v);
81  UPWIND3 interp_v_UPW3_v(v, upw_frac_v);
82  UPWIND3 interp_w_UPW3_v(w, upw_frac_v);
83 
84  bool u_is_regular = (u_type == amrex::FabType::regular);
85  bool v_is_regular = (v_type == amrex::FabType::regular);
86  bool w_is_regular = (w_type == amrex::FabType::regular);
87 
88  bool u_is_singlevalued = (u_type == amrex::FabType::singlevalued);
89  bool v_is_singlevalued = (v_type == amrex::FabType::singlevalued);
90  bool w_is_singlevalued = (w_type == amrex::FabType::singlevalued);
91 
92  // ---------------------------------------------------- //
93  // x-momentum
94  // ---------------------------------------------------- //
95 
96  amrex::ParallelFor(bxx_grown[0], bxx_grown[1], bxx_grown[2],
97  // x-momentum: x-face
98  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
99  {
100  bool not_covered = u_is_regular || (u_is_singlevalued && u_afrac_x(i,j,k) > zero);
101 
102  if ( not_covered ) {
103 
104  amrex::Real rho_u_avg = myhalf * (rho_u(i-1,j,k) * mf_uy_inv(i-1,j,0) + rho_u(i,j,k) * mf_uy_inv(i,j,0));
105 
106  int icell = 0;
107 
108  if (rho_u_avg>zero)
109  {
110  for (int ii=0; ii<ncells_u_horz; ++ii) {
111  if (u_cflag(i-ii-1,j,k).isCovered()) {
112  break;
113  }
114  icell++;
115  }
116  }
117  else if (rho_u_avg<zero)
118  {
119  for (int ii=0; ii<ncells_u_horz; ++ii) {
120  if (u_cflag(i+ii,j,k).isCovered()) {
121  break;
122  }
123  icell++;
124  }
125  }
126 
127  // Interpolate u using the highest-order scheme
128 
129  amrex::Real interpx = zero;
130  if (icell==ncells_u_horz) {
131  interp_u_horz.InterpolateInX(i,j,k,0,interpx,rho_u_avg);
132  } else {
133  if (icell<=1) {
134  interp_u_CEN2.InterpolateInX(i,j,k,0,interpx,rho_u_avg);
135  } else if (icell==2) {
136  interp_u_UPW3_horz.InterpolateInX(i,j,k,0,interpx,rho_u_avg);
137  }
138  }
139  flx_u_arr[0](i,j,k) = rho_u_avg * interpx;
140 
141  } else {
142 
143  flx_u_arr[0](i,j,k) = zero;
144  }
145  },
146  // x-momentum: y-face
147  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
148  {
149  bool not_covered = u_is_regular || (u_is_singlevalued && u_afrac_y(i,j,k) > zero);
150 
151  if ( not_covered ) {
152 
153  amrex::Real rho_v_avg = myhalf * (rho_v(i,j,k) * mf_vx_inv(i,j,0) + rho_v(i-1,j,k) * mf_vx_inv(i-1,j,0));
154 
155  int jcell = 0;
156 
157  if (rho_v_avg>zero)
158  {
159  for (int jj=0; jj<ncells_u_horz; ++jj) {
160  if (u_cflag(i,j-jj-1,k).isCovered()) {
161  break;
162  }
163  jcell++;
164  }
165  }
166  else if (rho_v_avg<zero)
167  {
168  for (int jj=0; jj<ncells_u_horz; ++jj) {
169  if (u_cflag(i,j+jj,k).isCovered()) {
170  break;
171  }
172  jcell++;
173  }
174  }
175 
176  // Interpolate u using the highest-order scheme
177 
178  amrex::Real interpy = zero;
179  if (jcell==ncells_u_horz) {
180  interp_u_horz.InterpolateInY(i,j,k,0,interpy,rho_v_avg);
181  } else {
182  if (jcell<=1) {
183  interp_u_CEN2.InterpolateInY(i,j,k,0,interpy,rho_v_avg);
184  } else if (jcell==2) {
185  interp_u_UPW3_horz.InterpolateInY(i,j,k,0,interpy,rho_v_avg);
186  }
187  }
188  flx_u_arr[1](i,j,k) = rho_v_avg * interpy;
189 
190  } else {
191 
192  flx_u_arr[1](i,j,k) = zero;
193  }
194  },
195  // x-momentum: z-face
196  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
197  {
198  bool not_covered = u_is_regular || (u_is_singlevalued && u_afrac_z(i,j,k) > zero);
199 
200  if ( not_covered ) {
201 
202  amrex::Real rho_w_avg = myhalf * (rho_w(i,j,k) + rho_w(i-1,j,k));
203 
204  int kcell = 0;
205 
206  if (rho_w_avg>zero)
207  {
208  for (int kk=0; kk<ncells_u_vert; ++kk) {
209  if (u_cflag(i,j,k-kk-1).isCovered()) {
210  break;
211  }
212  kcell++;
213  }
214  }
215  else if (rho_w_avg<zero)
216  {
217  for (int kk=0; kk<ncells_u_vert; ++kk) {
218  if (u_cflag(i,j,k+kk).isCovered()) {
219  break;
220  }
221  kcell++;
222  }
223  }
224 
225  // Interpolate u using the highest-order scheme
226 
227  amrex::Real interpz = zero;
228  if (kcell==ncells_u_vert) {
229  interp_u_vert.InterpolateInZ(i,j,k,0,interpz,rho_w_avg);
230  } else {
231  if (kcell<=1) {
232  interp_u_CEN2.InterpolateInZ(i,j,k,0,interpz,rho_w_avg);
233  } else if (kcell==2) {
234  interp_u_UPW3_vert.InterpolateInZ(i,j,k,0,interpz,rho_w_avg);
235  }
236  }
237  flx_u_arr[2](i,j,k) = rho_w_avg * interpz;
238 
239  } else {
240 
241  flx_u_arr[2](i,j,k) = zero;
242  }
243  });
244 
245  // ---------------------------------------------------- //
246  // y-momentum
247  // ---------------------------------------------------- //
248 
249  amrex::ParallelFor(bxy_grown[0], bxy_grown[1], bxy_grown[2],
250  // y-momentum: x-face
251  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
252  {
253  bool not_covered = v_is_regular || (v_is_singlevalued && v_afrac_x(i,j,k) > zero);
254 
255  if ( not_covered ) {
256 
257  amrex::Real rho_u_avg = myhalf * (rho_u(i,j,k) * mf_uy_inv(i,j,0) + rho_u(i,j-1,k) * mf_uy_inv(i,j-1,0));
258 
259  int icell = 0;
260 
261  if (rho_u_avg>zero)
262  {
263  for (int ii=0; ii<ncells_v_horz; ++ii) {
264  if (v_cflag(i-ii-1,j,k).isCovered()) {
265  break;
266  }
267  icell++;
268  }
269  }
270  else if (rho_u_avg<zero)
271  {
272  for (int ii=0; ii<ncells_v_horz; ++ii) {
273  if (v_cflag(i+ii,j,k).isCovered()) {
274  break;
275  }
276  icell++;
277  }
278  }
279 
280  // Interpolate v using the highest-order scheme
281 
282  amrex::Real interpx = zero;
283  if (icell==ncells_v_horz) {
284  interp_v_horz.InterpolateInX(i,j,k,0,interpx,rho_u_avg);
285  } else {
286  if (icell<=1) {
287  interp_v_CEN2.InterpolateInX(i,j,k,0,interpx,rho_u_avg);
288  } else if (icell==2) {
289  interp_v_UPW3_horz.InterpolateInX(i,j,k,0,interpx,rho_u_avg);
290  }
291  }
292  flx_v_arr[0](i,j,k) = rho_u_avg * interpx;
293 
294  } else {
295 
296  flx_v_arr[0](i,j,k) = zero;
297  }
298  },
299  // y-momentum: y-face
300  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
301  {
302  bool not_covered = v_is_regular || (v_is_singlevalued && v_afrac_y(i,j,k) > zero);
303 
304  if ( not_covered ) {
305 
306  amrex::Real rho_v_avg = myhalf * (rho_v(i,j,k) * mf_vx_inv(i,j,0) + rho_v(i,j-1,k) * mf_vx_inv(i,j-1,0));
307 
308  int jcell = 0;
309 
310  if (rho_v_avg>zero)
311  {
312  for (int jj=0; jj<ncells_v_horz; ++jj) {
313  if (v_cflag(i,j-jj-1,k).isCovered()) {
314  break;
315  }
316  jcell++;
317  }
318  }
319  else if (rho_v_avg<zero)
320  {
321  for (int jj=0; jj<ncells_v_horz; ++jj) {
322  if (v_cflag(i,j+jj,k).isCovered()) {
323  break;
324  }
325  jcell++;
326  }
327  }
328 
329  // Interpolate v using the highest-order scheme
330 
331  amrex::Real interpy = zero;
332  if (jcell==ncells_v_horz) {
333  interp_v_horz.InterpolateInY(i,j,k,0,interpy,rho_v_avg);
334  } else {
335  if (jcell<=1) {
336  interp_v_CEN2.InterpolateInY(i,j,k,0,interpy,rho_v_avg);
337  } else if (jcell==2) {
338  interp_v_UPW3_horz.InterpolateInY(i,j,k,0,interpy,rho_v_avg);
339  }
340  }
341  flx_v_arr[1](i,j,k) = rho_v_avg * interpy;
342 
343  } else {
344 
345  flx_v_arr[1](i,j,k) = zero;
346  }
347  },
348  // y-momentum: z-face
349  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
350  {
351  bool not_covered = v_is_regular || (v_is_singlevalued && v_afrac_z(i,j,k) > zero);
352 
353  if ( not_covered ) {
354 
355  amrex::Real rho_w_avg = myhalf * (rho_w(i,j,k) + rho_w(i,j-1,k));
356 
357  int kcell = 0;
358 
359  if (rho_w_avg>zero)
360  {
361  for (int kk=0; kk<ncells_v_vert; ++kk) {
362  if (v_cflag(i,j,k-kk-1).isCovered()) {
363  break;
364  }
365  kcell++;
366  }
367  }
368  else if (rho_w_avg<zero)
369  {
370  for (int kk=0; kk<ncells_v_vert; ++kk) {
371  if (v_cflag(i,j,k+kk).isCovered()) {
372  break;
373  }
374  kcell++;
375  }
376  }
377 
378  // Interpolate v using the highest-order scheme
379 
380  amrex::Real interpz = zero;
381  if (kcell==ncells_v_vert) {
382  interp_v_vert.InterpolateInZ(i,j,k,0,interpz,rho_w_avg);
383  } else {
384  if (kcell<=1) {
385  interp_v_CEN2.InterpolateInZ(i,j,k,0,interpz,rho_w_avg);
386  } else if (kcell==2) {
387  interp_v_UPW3_vert.InterpolateInZ(i,j,k,0,interpz,rho_w_avg);
388  }
389  }
390  flx_v_arr[2](i,j,k) = rho_w_avg * interpz;
391 
392  } else {
393 
394  flx_v_arr[2](i,j,k) = zero;
395  }
396  });
397 
398  // ---------------------------------------------------- //
399  // z-momentum
400  // ---------------------------------------------------- //
401 
402  amrex::ParallelFor(bxz_grown[0], bxz_grown[1], bxz_grown[2],
403  // z-momentum: x-face
404  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
405  {
406  bool not_covered = w_is_regular || (w_is_singlevalued && w_afrac_x(i,j,k) > zero);
407 
408  if ( not_covered ) {
409 
410  amrex::Real rho_u_avg = myhalf * (rho_u(i,j,k) + rho_u(i,j,k-1)) * mf_uy_inv(i,j,0);
411 
412  int icell = 0;
413 
414  if (rho_u_avg>zero)
415  {
416  for (int ii=0; ii<ncells_w_horz; ++ii) {
417  if (w_cflag(i-ii-1,j,k).isCovered()) {
418  break;
419  }
420  icell++;
421  }
422  }
423  else if (rho_u_avg<zero)
424  {
425  for (int ii=0; ii<ncells_w_horz; ++ii) {
426  if (w_cflag(i+ii,j,k).isCovered()) {
427  break;
428  }
429  icell++;
430  }
431  }
432 
433  // Interpolate w using the highest-order scheme
434 
435  amrex::Real interpx = zero;
436  if (icell==ncells_w_horz) {
437  interp_w_horz.InterpolateInX(i,j,k,0,interpx,rho_u_avg);
438  } else {
439  if (icell<=1) {
440  interp_w_CEN2.InterpolateInX(i,j,k,0,interpx,rho_u_avg);
441  } else if (icell==2) {
442  interp_w_UPW3_horz.InterpolateInX(i,j,k,0,interpx,rho_u_avg);
443  }
444  }
445  flx_w_arr[0](i,j,k) = rho_u_avg * interpx;
446 
447  } else {
448 
449  flx_w_arr[0](i,j,k) = zero;
450  }
451  },
452  // z-momentum: y-face
453  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
454  {
455  bool not_covered = w_is_regular || (w_is_singlevalued && w_afrac_y(i,j,k) > zero);
456 
457  if ( not_covered ) {
458 
459  amrex::Real rho_v_avg = myhalf * (rho_v(i,j,k) + rho_v(i,j,k-1)) * mf_vx_inv(i,j,0);
460 
461  int jcell = 0;
462 
463  if (rho_v_avg>zero)
464  {
465  for (int jj=0; jj<ncells_w_horz; ++jj) {
466  if (w_cflag(i,j-jj-1,k).isCovered()) {
467  break;
468  }
469  jcell++;
470  }
471  }
472  else if (rho_v_avg<zero)
473  {
474  for (int jj=0; jj<ncells_w_horz; ++jj) {
475  if (w_cflag(i,j+jj,k).isCovered()) {
476  break;
477  }
478  jcell++;
479  }
480  }
481 
482  // Interpolate w using the highest-order scheme
483 
484  amrex::Real interpy = zero;
485  if (jcell==ncells_w_horz) {
486  interp_w_horz.InterpolateInY(i,j,k,0,interpy,rho_v_avg);
487  } else {
488  if (jcell<=1) {
489  interp_w_CEN2.InterpolateInY(i,j,k,0,interpy,rho_v_avg);
490  } else if (jcell==2) {
491  interp_w_UPW3_horz.InterpolateInY(i,j,k,0,interpy,rho_v_avg);
492  }
493  }
494  flx_w_arr[1](i,j,k) = rho_v_avg * interpy;
495 
496  } else {
497 
498  flx_w_arr[1](i,j,k) = zero;
499  }
500  },
501  // z-momentum: z-face
502  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
503  {
504  bool not_covered = w_is_regular || (w_is_singlevalued && w_afrac_z(i,j,k) > zero);
505 
506  if ( not_covered ) {
507 
508  if (k==lo_z_face || k==hi_z_face+1) {
509 
510  flx_w_arr[2](i,j,k) = rho_w(i,j,k) * w(i,j,k);
511 
512  } else {
513  amrex::Real rho_w_avg = myhalf * (rho_w(i,j,k) + rho_w(i,j,k-1));
514 
515  int kcell = 0;
516 
517  if (rho_w_avg>zero)
518  {
519  for (int kk=0; kk<ncells_w_vert; ++kk) {
520  if (w_cflag(i,j,k-kk-1).isCovered()) {
521  break;
522  }
523  kcell++;
524  }
525  }
526  else if (rho_w_avg<zero)
527  {
528  for (int kk=0; kk<ncells_w_vert; ++kk) {
529  if (w_cflag(i,j,k+kk).isCovered()) {
530  break;
531  }
532  kcell++;
533  }
534  }
535 
536  // Adjust available cell numbers near the top and bottom boundaries.
537  // As in the non-EB operator, the second face from the wall uses
538  // 4th-order centered when the wall, not the EB, limits the stencil.
539 
540  bool wall_cen4 = false;
541  if (k == hi_z_face || k == lo_z_face+1) {
542  kcell=std::min(kcell,1);
543  } else if (k == hi_z_face-1 || k == lo_z_face+2 ) {
544  wall_cen4 = (kcell > 2);
545  kcell=std::min(kcell,2);
546  }
547 
548  // Interpolate w using the highest-order scheme
549 
550  amrex::Real interpz = zero;
551  if (kcell==ncells_w_vert) {
552  interp_w_vert.InterpolateInZ(i,j,k,0,interpz,rho_w_avg);
553  } else if (wall_cen4) {
554  interp_w_CEN4.InterpolateInZ(i,j,k,0,interpz,rho_w_avg);
555  } else {
556  if (kcell<=1) {
557  interp_w_CEN2.InterpolateInZ(i,j,k,0,interpz,rho_w_avg);
558  } else if (kcell==2) {
559  interp_w_UPW3_vert.InterpolateInZ(i,j,k,0,interpz,rho_w_avg);
560  }
561  }
562  flx_w_arr[2](i,j,k) = rho_w_avg * interpz;
563  } // lo_z_face, hi_z_face
564 
565  } else {
566 
567  flx_w_arr[2](i,j,k) = zero;
568  }
569  });
570 
571 }
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);})
zero
constexpr amrex::Real zero
Definition: ERF_NumericalConstants.H:36
myhalf
constexpr amrex::Real myhalf
Definition: ERF_NumericalConstants.H:41
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:19
CENTERED2
Definition: ERF_Interpolation_UPW.H:10
CENTERED4
Definition: ERF_Interpolation_UPW.H:335
UPWIND3
Definition: ERF_Interpolation_UPW.H:95
Here is the call graph for this function: